diff --git a/.clang-format b/.clang-format
index 24437d09..9215a1c4 100644
--- a/.clang-format
+++ b/.clang-format
@@ -78,13 +78,13 @@ BracedInitializerIndentWidth: 4
 
 IncludeBlocks: Regroup
 IncludeCategories:
-  - Regex: '^".*"'                # Internal includes inside double quotes
+  - Regex: '^".*"'                    # Internal includes inside double quotes
     Priority: 1
-  - Regex: '^<gl/[^>]+\.h(pp)?>'  # Internal includes inside angle brackets
+  - Regex: '^<(h)?gl/[^>]+\.h(pp)?>'  # Internal includes inside angle brackets
     Priority: 2
-  - Regex: '<.*\.h(pp)?>'         # External non-standard includes
+  - Regex: '<.*\.h(pp)?>'             # External non-standard includes
     Priority: 3
-  - Regex: '<.*>'                 # Standard library includes
+  - Regex: '<.*>'                     # Standard library includes
     Priority: 4
 SortIncludes:          CaseSensitive
 SortUsingDeclarations: LexicographicNumeric
diff --git a/.github/workflows/benchmarks.yaml b/.github/workflows/benchmarks.yaml
index af77b126..2f23bce9 100644
--- a/.github/workflows/benchmarks.yaml
+++ b/.github/workflows/benchmarks.yaml
@@ -42,4 +42,6 @@ jobs:
         run: |
           ./build_bench/benchmarks/cpp-gl-bench \
             --benchmark_repetitions=1 --benchmark_display_aggregates_only=true \
-            --bip-v 100
+            --bip-v 100 \
+            --hg-bfs-l --hg-bfs-al --hg-bfs-m --hg-bfs-e 100 --hg-bfs-esize 5 --hg-bfs-stride 2 \
+            --hg-b-bfs-l --hg-b-bfs-al --hg-b-bfs-m --hg-b-bfs-e 100 --hg-b-bfs-layer-width 5 --hg-b-bfs-stride 2
diff --git a/README.md b/README.md
index d431e53c..cff4e06b 100644
--- a/README.md
+++ b/README.md
@@ -160,7 +160,13 @@ python scripts/format.py -m -exe clang-format-18
 
 ### Building the Documentation
 
-The documentation is generated using Doxygen (for XML extraction), a custom Python script for concept parsing, and MkDocs (via `mike` for versioning). The process is automated via a Makefile.
+The documentation build process utilizes the following toolchain:
+
+- **[Doxygen](https://www.doxygen.nl/)**: Extracts the initial C++ API structure and documentation into XML format.
+- **Concept Parser**: A [custom Python script](/docs/scripts/gen_concept_docs.py) processes the XML to generate dedicated concept documentation pages.
+- **[MkDoxy (Custom Fork)](https://github.com/SpectraL519/MkDoxy)**: Injects the Doxygen XML data directly into the MkDocs build lifecycle.
+- **[MkDocs](https://www.mkdocs.org/) & [Material Theme](https://squidfunk.github.io/mkdocs-material/)**: Renders the final, searchable static HTML website.
+- **[mike](https://github.com/jimporter/mike)**: Manages versioning, allowing multiple versions of the documentation to coexist and be deployed simultaneously.
 
 To build and serve the documentation locally:
 
diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt
index 7ac49a7a..c82dbcf5 100644
--- a/benchmarks/CMakeLists.txt
+++ b/benchmarks/CMakeLists.txt
@@ -34,6 +34,8 @@ add_executable(cpp-gl-bench
     source/main.cpp
     source/runner.cpp
     suites/is_bipartite.cpp
+    suites/hg_bfs.cpp
+    suites/hg_b_bfs.cpp
 )
 
 target_compile_options(cpp-gl-bench PRIVATE
diff --git a/benchmarks/README.md b/benchmarks/README.md
index c6fdbfd9..4067292e 100644
--- a/benchmarks/README.md
+++ b/benchmarks/README.md
@@ -30,7 +30,7 @@ cmake -B build_bench -DBUILD_BENCHMARKS=ON -DBENCH_INCLUDE_BGL=ON -DCMAKE_BUILD_
 ```
 
 > [!NOTE]
-> 
+>
 > BGL comparative benchmarks are only defined for specific benchmark suites where a direct equivalent exists in Boost.
 
 <br/>
diff --git a/benchmarks/suites/hg_b_bfs.cpp b/benchmarks/suites/hg_b_bfs.cpp
new file mode 100644
index 00000000..fb4c969e
--- /dev/null
+++ b/benchmarks/suites/hg_b_bfs.cpp
@@ -0,0 +1,294 @@
+#include "gl/algorithm/core.hpp"
+#include "hgl/directional_tags.hpp"
+#include "runner.hpp"
+#include "suite.hpp"
+
+#include <gl/algorithm.hpp>
+#include <gl/graph.hpp>
+#include <hgl/algorithm.hpp>
+#include <hgl/conversion.hpp>
+#include <hgl/hypergraph.hpp>
+
+#include <benchmark/benchmark.h>
+#include <sys/types.h>
+
+namespace gl_bench::hg_b_bfs {
+
+// --- Hypergraph Topology Generator ---
+
+template <hgl::traits::c_bf_directed_hypergraph HypergraphType>
+HypergraphType gen_bf_overlapping_chain_hypergraph(
+    const std::size_t n_hyperedges, const std::size_t layer_width, const std::size_t stride
+) {
+    using id_type = typename HypergraphType::id_type;
+
+    // V = maximum index reached by the last hyperedge's head.
+    const auto n_vertices = static_cast<id_type>((n_hyperedges - 1) * stride + 2 * layer_width);
+
+    HypergraphType hgraph{n_vertices, static_cast<id_type>(n_hyperedges)};
+
+    for (id_type e = 0; e < static_cast<id_type>(n_hyperedges); ++e) {
+        const auto tail_start = static_cast<id_type>(e * stride);
+        const auto tail_end = tail_start + static_cast<id_type>(layer_width);
+
+        const auto head_start = tail_end;
+        const auto head_end = head_start + static_cast<id_type>(layer_width);
+
+        hgraph.bind_tail(std::views::iota(tail_start, tail_end), e);
+        hgraph.bind_head(std::views::iota(head_start, head_end), e);
+    }
+
+    return hgraph;
+}
+
+template <hgl::traits::c_bf_directed_hypergraph HypergraphType>
+requires hgl::traits::c_flat_list_hypergraph<HypergraphType>
+HypergraphType gen_bf_overlapping_chain_hypergraph(
+    const std::size_t n_hyperedges, const std::size_t layer_width, const std::size_t stride
+) {
+    using traits_type = typename HypergraphType::traits_type;
+    using layout_tag = typename traits_type::layout_tag;
+
+    using list_repr_tag = hgl::repr::list_t<layout_tag>;
+    using list_hypergraph = hgl::traits::swap_repr_tag_t<HypergraphType, list_repr_tag>;
+
+    return hgl::to<typename traits_type::representation_tag>(
+        gen_bf_overlapping_chain_hypergraph<list_hypergraph>(n_hyperedges, layer_width, stride)
+    );
+}
+
+// --- Native HGL Backward BFS Benchmark ---
+
+template <hgl::traits::c_bf_directed_hypergraph Hypergraph>
+void bm_hgl_backward_bfs(benchmark::State& state) {
+    using id_type = typename Hypergraph::id_type;
+
+    const auto n_hedges = static_cast<std::size_t>(state.range(0));
+    const auto layer_width = static_cast<std::size_t>(state.range(1));
+    const auto stride = static_cast<std::size_t>(state.range(2));
+
+    auto hg = gen_bf_overlapping_chain_hypergraph<Hypergraph>(n_hedges, layer_width, stride);
+
+    auto roots = std::views::iota(id_type{0}, static_cast<id_type>(layer_width))
+               | std::ranges::to<std::vector<id_type>>();
+
+    for (auto _ : state) {
+        auto search_tree = hgl::algorithm::backward_bfs(hg, roots);
+        benchmark::DoNotOptimize(search_tree);
+    }
+
+    state.counters["Vertices"] = static_cast<double>(hg.n_vertices());
+    state.counters["Hyperedges"] = static_cast<double>(hg.n_hyperedges());
+    state.counters["Incidences"] = static_cast<double>(n_hedges * layer_width * 2);
+}
+
+// --- Incidence Graph Utilities ---
+
+/// @brief Executes a Backward BFS equivalent on an Incidence Graph.
+template <gl::traits::c_directed_graph IncidenceGraph>
+bool incidence_backward_bfs(
+    const IncidenceGraph& ig,
+    const std::vector<typename IncidenceGraph::id_type>& roots,
+    const typename IncidenceGraph::id_type original_n_vertices
+) {
+    using id_type = typename IncidenceGraph::id_type;
+
+    std::vector<bool> visited_v(original_n_vertices, false);
+    auto tail_unvisited =
+        ig.in_degree_map() | std::views::drop(original_n_vertices)
+        | std::ranges::to<std::vector<gl::size_type>>();
+
+    auto root_nodes =
+        roots | std::views::transform([](const id_type root_id) {
+            return gl::algorithm::search_node<IncidenceGraph>{root_id};
+        })
+        | std::ranges::to<std::vector>();
+
+    auto visit_vertex_pred = [&](id_type v) {
+        if (v < original_n_vertices)
+            return not visited_v[gl::to_idx(v)];
+        return true;
+    };
+
+    auto visit = [&](id_type v, id_type /*p*/) {
+        if (v < original_n_vertices)
+            visited_v[gl::to_idx(v)] = true;
+        return true;
+    };
+
+    auto enqueue_node_pred =
+        [&](id_type target_id, const auto& /*edge*/) -> gl::algorithm::decision {
+        if (target_id >= original_n_vertices) {
+            const auto he_idx = target_id - original_n_vertices;
+            return --tail_unvisited[gl::to_idx(he_idx)] == 0uz;
+        }
+        else {
+            return not visited_v[gl::to_idx(target_id)];
+        }
+    };
+
+    return gl::algorithm::bfs(ig, root_nodes, visit_vertex_pred, visit, enqueue_node_pred);
+}
+
+// --- GL Incidence Graph Backward BFS Benchmark ---
+
+template <
+    hgl::traits::c_bf_directed_hypergraph Hypergraph,
+    gl::traits::c_directed_graph IncidenceGraph>
+void bm_gl_incidence_backward_bfs(benchmark::State& state) {
+    using id_type = typename IncidenceGraph::id_type;
+
+    const auto n_hedges = static_cast<std::size_t>(state.range(0));
+    const auto layer_width = static_cast<std::size_t>(state.range(1));
+    const auto stride = static_cast<std::size_t>(state.range(2));
+
+    if constexpr (gl::traits::c_adjacency_matrix_graph<IncidenceGraph>) {
+        const auto n_vertices = (n_hedges * stride) + layer_width;
+        const auto ig_vertices = n_vertices + n_hedges;
+
+        // Hardcoded safety limit for ~16 GB of RAM (V_ig = 65,000)
+        if (ig_vertices > 65000) {
+            state.SkipWithError("Matrix requires > 16.0 GB of memory; skipping.");
+            return;
+        }
+    }
+
+    auto hg = gen_bf_overlapping_chain_hypergraph<Hypergraph>(n_hedges, layer_width, stride);
+    auto ig = hgl::incidence_graph<IncidenceGraph>(hg);
+
+    auto roots = std::views::iota(id_type{0}, static_cast<id_type>(layer_width))
+               | std::ranges::to<std::vector<id_type>>();
+
+    for (auto _ : state) {
+        bool completed = incidence_backward_bfs(ig, roots, static_cast<id_type>(hg.n_vertices()));
+        benchmark::DoNotOptimize(completed);
+    }
+
+    state.counters["IG_Vertices"] = static_cast<double>(ig.n_vertices());
+    state.counters["IG_Edges"] = static_cast<double>(ig.n_edges());
+}
+
+// --- Suite Setup & Registration ---
+
+void add_args(argon::argument_parser& parser) {
+    auto& group = parser.add_group("Backward BFS Benchmark Suite Options (hg-b-bfs)");
+    parser.add_optional_argument<std::size_t>(group, "hg-b-bfs-e")
+        .default_values(1000uz)
+        .help("Number of BF-directed hyperedges");
+    parser.add_optional_argument<std::size_t>(group, "hg-b-bfs-layer-width")
+        .default_values(10uz)
+        .help("Number of vertices per layer (tail and head sizes)");
+    parser.add_optional_argument<std::size_t>(group, "hg-b-bfs-stride")
+        .default_values(2uz)
+        .help("Vertex shift between consecutive hyperedges (smaller = denser)");
+    parser.add_flag("hg-b-bfs-l").help("Execute the benchmark for list models (bidirectional)");
+    parser.add_flag("hg-b-bfs-al")
+        .help("Execute the benchmark for asymmetric list models (WARNING: High execution times "
+              "expected)");
+    parser.add_flag("hg-b-bfs-m").help("Execute the benchmark for matrix models");
+}
+
+void register_benchmarks(const argon::argument_parser& parser) {
+    const auto n_hedges = static_cast<int64_t>(parser.value<std::size_t>("hg-b-bfs-e"));
+    const auto layer_width =
+        static_cast<int64_t>(parser.value<std::size_t>("hg-b-bfs-layer-width"));
+    const auto stride = static_cast<int64_t>(parser.value<std::size_t>("hg-b-bfs-stride"));
+
+    using gl_list = gl::list_graph<gl::directed_t>;
+    using gl_flat_list = gl::flat_list_graph<gl::directed_t>;
+    using gl_matrix = gl::matrix_graph<gl::directed_t>;
+    using gl_flat_matrix = gl::flat_matrix_graph<gl::directed_t>;
+
+    using hgl_list = hgl::list_hypergraph<hgl::repr::bidirectional_t, hgl::bf_directed_t>;
+    using hgl_flat_list = hgl::flat_list_hypergraph<hgl::repr::bidirectional_t, hgl::bf_directed_t>;
+
+    using hgl_v_list = hgl::list_hypergraph<hgl::repr::vertex_major_t, hgl::bf_directed_t>;
+    using hgl_v_flat_list =
+        hgl::flat_list_hypergraph<hgl::repr::vertex_major_t, hgl::bf_directed_t>;
+
+    using hgl_e_list = hgl::list_hypergraph<hgl::repr::hyperedge_major_t, hgl::bf_directed_t>;
+    using hgl_e_flat_list =
+        hgl::flat_list_hypergraph<hgl::repr::hyperedge_major_t, hgl::bf_directed_t>;
+
+    using hgl_v_matrix = hgl::matrix_hypergraph<hgl::repr::vertex_major_t, hgl::bf_directed_t>;
+    using hgl_e_matrix = hgl::matrix_hypergraph<hgl::repr::hyperedge_major_t, hgl::bf_directed_t>;
+
+    using hgl_v_flat_matrix =
+        hgl::flat_matrix_hypergraph<hgl::repr::vertex_major_t, hgl::bf_directed_t>;
+    using hgl_e_flat_matrix =
+        hgl::flat_matrix_hypergraph<hgl::repr::hyperedge_major_t, hgl::bf_directed_t>;
+
+    if (parser.value<bool>("hg-b-bfs-l")) {
+        benchmark::RegisterBenchmark("b_bfs/HGL/list/bidir", bm_hgl_backward_bfs<hgl_list>)
+            ->Args({n_hedges, layer_width, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("b_bfs/HGL/flat_list/bidir", bm_hgl_backward_bfs<hgl_flat_list>)
+            ->Args({n_hedges, layer_width, stride})
+            ->Unit(benchmark::kMillisecond);
+
+        benchmark::
+            RegisterBenchmark("b_bfs/INCIDENCE/list", bm_gl_incidence_backward_bfs<hgl_list, gl_list>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("b_bfs/INCIDENCE/flat_list", bm_gl_incidence_backward_bfs<hgl_list, gl_flat_list>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+    }
+
+    if (parser.value<bool>("hg-b-bfs-al")) {
+        benchmark::RegisterBenchmark("bfs/HGL/list/v_major", bm_hgl_backward_bfs<hgl_v_list>)
+            ->Args({n_hedges, layer_width, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("bfs/HGL/flat_list/v_major", bm_hgl_backward_bfs<hgl_v_flat_list>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+
+        benchmark::RegisterBenchmark("bfs/HGL/list/e_major", bm_hgl_backward_bfs<hgl_e_list>)
+            ->Args({n_hedges, layer_width, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("bfs/HGL/flat_list/e_major", bm_hgl_backward_bfs<hgl_e_flat_list>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+    }
+
+    if (parser.value<bool>("hg-b-bfs-m")) {
+        benchmark::RegisterBenchmark("b_bfs/HGL/matrix/v_major", bm_hgl_backward_bfs<hgl_v_matrix>)
+            ->Args({n_hedges, layer_width, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("b_bfs/HGL/matrix/e_major", bm_hgl_backward_bfs<hgl_e_matrix>)
+            ->Args({n_hedges, layer_width, stride})
+            ->Unit(benchmark::kMillisecond);
+
+        benchmark::
+            RegisterBenchmark("b_bfs/HGL/flat_matrix/v_major", bm_hgl_backward_bfs<hgl_v_flat_matrix>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("b_bfs/HGL/flat_matrix/e_major", bm_hgl_backward_bfs<hgl_e_flat_matrix>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+
+        benchmark::
+            RegisterBenchmark("b_bfs/INCIDENCE/matrix", bm_gl_incidence_backward_bfs<hgl_list, gl_matrix>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("b_bfs/INCIDENCE/flat_matrix", bm_gl_incidence_backward_bfs<hgl_list, gl_flat_matrix>)
+                ->Args({n_hedges, layer_width, stride})
+                ->Unit(benchmark::kMillisecond);
+    }
+}
+
+namespace {
+bool _registered = []() {
+    gl_bench::runner::get().add_suite(
+        "hg-b-bfs", suite{.add_args = add_args, .register_benchmarks = register_benchmarks}
+    );
+    return true;
+}();
+} // namespace
+
+} // namespace gl_bench::hg_b_bfs
diff --git a/benchmarks/suites/hg_bfs.cpp b/benchmarks/suites/hg_bfs.cpp
new file mode 100644
index 00000000..e0c76bd2
--- /dev/null
+++ b/benchmarks/suites/hg_bfs.cpp
@@ -0,0 +1,219 @@
+#include "hgl/repr/layout_tags.hpp"
+#include "runner.hpp"
+#include "suite.hpp"
+
+#include <gl/algorithm.hpp>
+#include <gl/graph.hpp>
+#include <hgl/algorithm.hpp>
+#include <hgl/conversion.hpp>
+#include <hgl/hypergraph.hpp>
+
+#include <benchmark/benchmark.h>
+
+namespace gl_bench::hg_bfs {
+
+// --- Hypergraph Topology Generator ---
+
+template <hgl::traits::c_undirected_hypergraph HypergraphType>
+HypergraphType gen_sliding_window_hypergraph(
+    const std::size_t n_hyperedges, const std::size_t degree, const std::size_t stride
+) {
+    using id_type = typename HypergraphType::id_type;
+
+    // V = (E - 1) * stride + degree
+    const auto n_vertices = static_cast<id_type>((n_hyperedges - 1) * stride + degree);
+    HypergraphType hgraph{n_vertices, static_cast<id_type>(n_hyperedges)};
+
+    for (id_type e = 0; e < static_cast<id_type>(n_hyperedges); ++e) {
+        const auto start_v = static_cast<id_type>(e * stride);
+        const auto end_v = start_v + static_cast<id_type>(degree);
+        hgraph.bind(std::views::iota(start_v, end_v), e);
+    }
+
+    return hgraph;
+}
+
+template <hgl::traits::c_undirected_hypergraph HypergraphType>
+requires hgl::traits::c_flat_list_hypergraph<HypergraphType>
+HypergraphType gen_sliding_window_hypergraph(
+    const std::size_t n_hyperedges, const std::size_t degree, const std::size_t stride
+) {
+    using traits_type = typename HypergraphType::traits_type;
+    using layout_tag = typename traits_type::layout_tag;
+
+    using list_repr_tag = hgl::repr::list_t<layout_tag>;
+    using list_hypergraph = hgl::traits::swap_repr_tag_t<HypergraphType, list_repr_tag>;
+
+    return hgl::to<typename traits_type::representation_tag>(
+        gen_sliding_window_hypergraph<list_hypergraph>(n_hyperedges, degree, stride)
+    );
+}
+
+// --- HGL Benchmark ---
+
+template <hgl::traits::c_undirected_hypergraph Hypergraph>
+void bm_hgl_bfs(benchmark::State& state) {
+    const auto n_hedges = static_cast<std::size_t>(state.range(0));
+    const auto he_size = static_cast<std::size_t>(state.range(1));
+    const auto stride = static_cast<std::size_t>(state.range(2));
+
+    auto hg = gen_sliding_window_hypergraph<Hypergraph>(n_hedges, he_size, stride);
+
+    for (auto _ : state) {
+        auto search_tree = hgl::algorithm::breadth_first_search(hg, 0uz);
+        benchmark::DoNotOptimize(search_tree);
+    }
+
+    state.counters["Vertices"] = static_cast<double>(hg.n_vertices());
+    state.counters["Hyperedges"] = static_cast<double>(hg.n_hyperedges());
+    state.counters["Incidences"] = static_cast<double>(n_hedges * he_size);
+}
+
+// --- GL Incidence Graph Benchmark ---
+
+template <
+    hgl::traits::c_undirected_hypergraph Hypergraph,
+    gl::traits::c_undirected_graph IncidenceGraph>
+void bm_gl_incidence_bfs(benchmark::State& state) {
+    const auto n_hedges = static_cast<std::size_t>(state.range(0));
+    const auto he_size = static_cast<std::size_t>(state.range(1));
+    const auto stride = static_cast<std::size_t>(state.range(2));
+
+    if constexpr (hgl::traits::c_incidence_matrix_hypergraph<Hypergraph>) {
+        const auto n_vertices = (n_hedges - 1) * stride + he_size;
+        const auto ig_vertices = n_vertices + n_hedges;
+
+        // Hardcoded safety limit for ~16 GB of RAM (V_ig = 65,000)
+        if (ig_vertices > 65000) {
+            state.SkipWithError("Matrix requires > 16.0 GB of memory; skipping.");
+            return;
+        }
+    }
+
+    auto hg = gen_sliding_window_hypergraph<Hypergraph>(n_hedges, he_size, stride);
+    auto ig = hgl::incidence_graph<IncidenceGraph>(hg);
+
+    for (auto _ : state) {
+        auto pred_map = gl::algorithm::breadth_first_search(ig, 0uz);
+        benchmark::DoNotOptimize(pred_map);
+    }
+
+    state.counters["IG_Vertices"] = static_cast<double>(ig.n_vertices());
+    state.counters["IG_Edges"] = static_cast<double>(ig.n_edges());
+}
+
+void add_args(argon::argument_parser& parser) {
+    auto& group = parser.add_group("Hypergraph BFS Benchmark Suite Options (hg-bfs)");
+    parser.add_optional_argument<std::size_t>(group, "hg-bfs-e")
+        .default_values(1000uz)
+        .help("Number of hyperedges");
+    parser.add_optional_argument<std::size_t>(group, "hg-bfs-esize")
+        .default_values(100uz)
+        .help("Hyperedge size");
+    parser.add_optional_argument<std::size_t>(group, "hg-bfs-stride")
+        .default_values(2uz)
+        .help("Vertex shift between consecutive hyperedges (smaller = denser)");
+    parser.add_flag("hg-bfs-l").help("Execute the benchmark for list models (bidirectional)");
+    parser.add_flag("hg-bfs-al")
+        .help("Execute the benchmark for asymmetric list models (WARNING: High execution times "
+              "expected)");
+    parser.add_flag("hg-bfs-m").help("Execute the benchmark for matrix models");
+}
+
+void register_benchmarks(const argon::argument_parser& parser) {
+    const auto n_hedges = static_cast<int64_t>(parser.value<std::size_t>("hg-bfs-e"));
+    const auto he_size = static_cast<int64_t>(parser.value<std::size_t>("hg-bfs-esize"));
+    const auto stride = static_cast<int64_t>(parser.value<std::size_t>("hg-bfs-stride"));
+
+    using gl_list = gl::list_graph<gl::undirected_t>;
+    using gl_flat_list = gl::flat_list_graph<gl::undirected_t>;
+    using gl_matrix = gl::matrix_graph<gl::undirected_t>;
+    using gl_flat_matrix = gl::flat_matrix_graph<gl::undirected_t>;
+
+    using hgl_list = hgl::list_hypergraph<hgl::repr::bidirectional_t, hgl::undirected_t>;
+    using hgl_flat_list = hgl::flat_list_hypergraph<hgl::repr::bidirectional_t, hgl::undirected_t>;
+
+    using hgl_v_list = hgl::list_hypergraph<hgl::repr::vertex_major_t, hgl::undirected_t>;
+    using hgl_v_flat_list = hgl::flat_list_hypergraph<hgl::repr::vertex_major_t, hgl::undirected_t>;
+
+    using hgl_e_list = hgl::list_hypergraph<hgl::repr::hyperedge_major_t, hgl::undirected_t>;
+    using hgl_e_flat_list =
+        hgl::flat_list_hypergraph<hgl::repr::hyperedge_major_t, hgl::undirected_t>;
+
+    using hgl_v_matrix = hgl::matrix_hypergraph<hgl::repr::vertex_major_t, hgl::undirected_t>;
+    using hgl_e_matrix = hgl::matrix_hypergraph<hgl::repr::hyperedge_major_t, hgl::undirected_t>;
+
+    using hgl_v_flat_matrix =
+        hgl::flat_matrix_hypergraph<hgl::repr::vertex_major_t, hgl::undirected_t>;
+    using hgl_e_flat_matrix =
+        hgl::flat_matrix_hypergraph<hgl::repr::hyperedge_major_t, hgl::undirected_t>;
+
+    if (parser.value<bool>("hg-bfs-l")) {
+        benchmark::RegisterBenchmark("bfs/HGL/list/bidir", bm_hgl_bfs<hgl_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("bfs/HGL/flat_list/bidir", bm_hgl_bfs<hgl_flat_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+
+        benchmark::RegisterBenchmark("bfs/INCIDENCE/list", bm_gl_incidence_bfs<hgl_list, gl_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("bfs/INCIDENCE/flat_list", bm_gl_incidence_bfs<hgl_list, gl_flat_list>)
+                ->Args({n_hedges, he_size, stride})
+                ->Unit(benchmark::kMillisecond);
+    }
+
+    if (parser.value<bool>("hg-bfs-al")) {
+        benchmark::RegisterBenchmark("bfs/HGL/list/v_major", bm_hgl_bfs<hgl_v_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("bfs/HGL/flat_list/v_major", bm_hgl_bfs<hgl_v_flat_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+
+        benchmark::RegisterBenchmark("bfs/HGL/list/e_major", bm_hgl_bfs<hgl_e_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("bfs/HGL/flat_list/e_major", bm_hgl_bfs<hgl_e_flat_list>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+    }
+
+    if (parser.value<bool>("hg-bfs-m")) {
+        benchmark::RegisterBenchmark("bfs/HGL/matrix/v_major", bm_hgl_bfs<hgl_v_matrix>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("bfs/HGL/matrix/e_major", bm_hgl_bfs<hgl_e_matrix>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+
+        benchmark::RegisterBenchmark("bfs/HGL/flat_matrix/v_major", bm_hgl_bfs<hgl_v_flat_matrix>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+        benchmark::RegisterBenchmark("bfs/HGL/flat_matrix/e_major", bm_hgl_bfs<hgl_e_flat_matrix>)
+            ->Args({n_hedges, he_size, stride})
+            ->Unit(benchmark::kMillisecond);
+
+        benchmark::
+            RegisterBenchmark("bfs/INCIDENCE/matrix", bm_gl_incidence_bfs<hgl_list, gl_matrix>)
+                ->Args({n_hedges, he_size, stride})
+                ->Unit(benchmark::kMillisecond);
+        benchmark::
+            RegisterBenchmark("bfs/INCIDENCE/flat_matrix", bm_gl_incidence_bfs<hgl_list, gl_flat_matrix>)
+                ->Args({n_hedges, he_size, stride})
+                ->Unit(benchmark::kMillisecond);
+    }
+}
+
+namespace {
+bool _registered = []() {
+    gl_bench::runner::get().add_suite(
+        "hg-bfs", suite{.add_args = add_args, .register_benchmarks = register_benchmarks}
+    );
+    return true;
+}();
+} // namespace
+
+} // namespace gl_bench::hg_bfs
diff --git a/benchmarks/suites/is_bipartite.cpp b/benchmarks/suites/is_bipartite.cpp
index c1857cc5..e637f85d 100644
--- a/benchmarks/suites/is_bipartite.cpp
+++ b/benchmarks/suites/is_bipartite.cpp
@@ -78,16 +78,10 @@ void register_benchmarks(const argon::argument_parser& parser) {
     const auto n_vertices = static_cast<int64_t>(parser.value<std::size_t>("bip-v"));
 
     // CPP-GL Adjacency List Benchmarks
-    using gl_list_u32 = gl::graph<gl::list_graph_traits<
-        gl::undirected_t,
-        gl::empty_properties,
-        gl::empty_properties,
-        std::uint32_t>>;
-    using gl_list_u64 = gl::graph<gl::list_graph_traits<
-        gl::undirected_t,
-        gl::empty_properties,
-        gl::empty_properties,
-        std::uint64_t>>;
+    using gl_list_u32 =
+        gl::list_graph<gl::undirected_t, gl::empty_properties, gl::empty_properties, std::uint32_t>;
+    using gl_list_u64 =
+        gl::list_graph<gl::undirected_t, gl::empty_properties, gl::empty_properties, std::uint64_t>;
 
     benchmark::RegisterBenchmark("is_bipartite/CPP-GL/list/u32", bm_gl_is_bipartite<gl_list_u32>)
         ->Arg(n_vertices)
@@ -97,16 +91,10 @@ void register_benchmarks(const argon::argument_parser& parser) {
         ->Unit(benchmark::kMillisecond);
 
     // CPP-GL Flat Adjacency List Benchmarks
-    using gl_flat_list_u32 = gl::graph<gl::flat_list_graph_traits<
-        gl::undirected_t,
-        gl::empty_properties,
-        gl::empty_properties,
-        std::uint32_t>>;
-    using gl_flat_list_u64 = gl::graph<gl::flat_list_graph_traits<
-        gl::undirected_t,
-        gl::empty_properties,
-        gl::empty_properties,
-        std::uint64_t>>;
+    using gl_flat_list_u32 = gl::
+        flat_list_graph<gl::undirected_t, gl::empty_properties, gl::empty_properties, std::uint32_t>;
+    using gl_flat_list_u64 = gl::
+        flat_list_graph<gl::undirected_t, gl::empty_properties, gl::empty_properties, std::uint64_t>;
 
     benchmark::
         RegisterBenchmark("is_bipartite/CPP-GL/flat_list/u32", bm_gl_is_bipartite<gl_flat_list_u32>)
@@ -118,16 +106,10 @@ void register_benchmarks(const argon::argument_parser& parser) {
             ->Unit(benchmark::kMillisecond);
 
     // CPP-GL Adjacency Matrix Benchmarks
-    using gl_matrix_u32 = gl::graph<gl::matrix_graph_traits<
-        gl::undirected_t,
-        gl::empty_properties,
-        gl::empty_properties,
-        std::uint32_t>>;
-    using gl_matrix_u64 = gl::graph<gl::matrix_graph_traits<
-        gl::undirected_t,
-        gl::empty_properties,
-        gl::empty_properties,
-        std::uint64_t>>;
+    using gl_matrix_u32 = gl::
+        matrix_graph<gl::undirected_t, gl::empty_properties, gl::empty_properties, std::uint32_t>;
+    using gl_matrix_u64 = gl::
+        matrix_graph<gl::undirected_t, gl::empty_properties, gl::empty_properties, std::uint64_t>;
 
     benchmark::RegisterBenchmark("is_bipartite/CPP-GL/matrix/u32", bm_gl_is_bipartite<gl_matrix_u32>)
         ->Arg(n_vertices)
@@ -137,16 +119,16 @@ void register_benchmarks(const argon::argument_parser& parser) {
         ->Unit(benchmark::kMillisecond);
 
     // CPP-GL Flat Adjacency Matrix Benchmarks
-    using gl_flat_matrix_u32 = gl::graph<gl::flat_matrix_graph_traits<
+    using gl_flat_matrix_u32 = gl::flat_matrix_graph<
         gl::undirected_t,
         gl::empty_properties,
         gl::empty_properties,
-        std::uint32_t>>;
-    using gl_flat_matrix_u64 = gl::graph<gl::flat_matrix_graph_traits<
+        std::uint32_t>;
+    using gl_flat_matrix_u64 = gl::flat_matrix_graph<
         gl::undirected_t,
         gl::empty_properties,
         gl::empty_properties,
-        std::uint64_t>>;
+        std::uint64_t>;
 
     benchmark::
         RegisterBenchmark("is_bipartite/CPP-GL/flat_matrix/u32", bm_gl_is_bipartite<gl_flat_matrix_u32>)
diff --git a/docs/gl/architecture.md b/docs/gl/architecture.md
index ac525521..8d63477a 100644
--- a/docs/gl/architecture.md
+++ b/docs/gl/architecture.md
@@ -48,7 +48,7 @@ Formally, a graph $G = (V, E)$ consists of a set of vertices $V$ and a set of ed
 
 - [**directed_t**](../cpp-gl/structgl_1_1directed__t.md) : Specifies a **directed graph** configuration where edges are defined as ordered pairs $(u, v)$ such that $u, v \in V$. In this graph type, a connection from vertex $u$ to vertex $v$ is structurally distinct from a connection from $v$ to $u$, and the existence of one directed edge does not imply the existence of the other.
 
-- [**gl::undirected_t**](../cpp-gl/structgl_1_1undirected__t.md) : Edges are unordered pairs $\{u, v\}$ where $u, v \in V$. The library automatically manages the bidirectional nature of these connections, ensuring that an edge between $u$ and $v$ is recognized during both traversal and structural degree calculations regardless of the order of endpoints.
+- [**undirected_t**](../cpp-gl/structgl_1_1undirected__t.md) : Edges are unordered pairs $\{u, v\}$ where $u, v \in V$. The library automatically manages the bidirectional nature of these connections, ensuring that an edge between $u$ and $v$ is recognized during both traversal and structural degree calculations regardless of the order of endpoints.
 
 ### IDs vs. Descriptors
 
diff --git a/include/gl/graph.hpp b/include/gl/graph.hpp
index a41a3e7d..d8830d7b 100644
--- a/include/gl/graph.hpp
+++ b/include/gl/graph.hpp
@@ -1143,7 +1143,7 @@ class graph final {
     /// @param is The source input stream.
     /// @param g The graph instance to populate.
     /// @return The stream reference for chaining.
-    friend gl_attr_force_inline std::istream& operator>>(std::istream& is, graph& g) {
+    gl_attr_force_inline friend std::istream& operator>>(std::istream& is, graph& g) {
         return g._gsf_read(is);
     }
 
diff --git a/include/gl/util/ranges.hpp b/include/gl/util/ranges.hpp
index 9fdd4b50..bcfa4dfd 100644
--- a/include/gl/util/ranges.hpp
+++ b/include/gl/util/ranges.hpp
@@ -65,7 +65,7 @@ template <std::ranges::forward_range R>
 ///
 /// > [!INFO] Time Complexity
 /// >
-/// > $O(N)$ where $N$ is the number of elements in the range.
+/// > \f$O(N)\f$ where $N$ is the number of elements in the range.
 template <std::ranges::forward_range R>
 [[nodiscard]] constexpr bool all_equal(
     R&& range, const std::ranges::range_value_t<R>& value
diff --git a/include/hgl/hypergraph.hpp b/include/hgl/hypergraph.hpp
index b18a57f7..0d9b89d2 100644
--- a/include/hgl/hypergraph.hpp
+++ b/include/hgl/hypergraph.hpp
@@ -144,7 +144,7 @@ struct to_impl;
 /// - **Standard Range Support**: Exposes lightweight views compliant with C++20 `std::ranges`, enabling functional-style iteration and algorithms.
 ///
 /// ### Basic Definitions
-/// A hypergraph \f$G = (V, E)\f$ consists of a set of vertices \f$V\f$ and a set of hyperedges \f$E\f$.
+/// A hypergraph \f$H = (V, E)\f$ consists of a set of vertices \f$V\f$ and a set of hyperedges \f$E\f$.
 ///
 /// - For undirected graphs, a hyperedge is a subset of the vertex set. Formally \f$E \subseteq 2^V\f$ and \f$e \in E \implies e \subseteq V\f$.
 /// - For BF-directed graphs, a hyperedge is an ordered pair of disjoint subsets of the vertex set - the *tail* (sources) and *head* (targets) of the hyperedge.
@@ -1975,7 +1975,7 @@ class hypergraph final {
     /// @param is The source input stream.
     /// @param g The hypergraph instance to populate.
     /// @return The stream reference for chaining.
-    friend gl_attr_force_inline std::istream& operator>>(std::istream& is, hypergraph& hg) {
+    gl_attr_force_inline friend std::istream& operator>>(std::istream& is, hypergraph& hg) {
         return hg._hgsf_read(is);
     }
 
diff --git a/include/hgl/impl/flat_incidence_list.hpp b/include/hgl/impl/flat_incidence_list.hpp
index e2e1c738..8a6c55c1 100644
--- a/include/hgl/impl/flat_incidence_list.hpp
+++ b/include/hgl/impl/flat_incidence_list.hpp
@@ -407,23 +407,37 @@ class flat_incidence_list<hgl::bf_directed_t, ReprTag> final {
 
     // --- binding methods ---
 
-    gl_attr_force_inline void bind_tail(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_tail(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
-        this->_remove_no_align(this->_head_storage, major_id, minor_id);
+
+        if (detail::contains(this->_head_storage[to_idx(major_id)], minor_id)) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the head of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         detail::unique_insert(this->_tail_storage, major_id, minor_id);
     }
 
-    gl_attr_force_inline void bind_head(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_head(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
-        this->_remove_no_align(this->_tail_storage, major_id, minor_id);
+
+        if (detail::contains(this->_tail_storage[to_idx(major_id)], minor_id)) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the tail of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         detail::unique_insert(this->_head_storage, major_id, minor_id);
     }
 
-    gl_attr_force_inline void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
+    void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
         this->_remove_no_align(this->_tail_storage, major_id, minor_id);
         this->_remove_no_align(this->_head_storage, major_id, minor_id);
@@ -807,14 +821,14 @@ class flat_incidence_list<DirectionalTag, ReprTag> final {
             return this->_e_list.are_bound(vertex_id, hyperedge_id);
     }
 
-    gl_attr_force_inline void bind_tail(const id_type vertex_id, const id_type hyperedge_id) noexcept
+    gl_attr_force_inline void bind_tail(const id_type vertex_id, const id_type hyperedge_id)
     requires std::same_as<DirectionalTag, hgl::bf_directed_t>
     {
         this->_v_list.bind_tail(vertex_id, hyperedge_id);
         this->_e_list.bind_tail(vertex_id, hyperedge_id);
     }
 
-    gl_attr_force_inline void bind_head(const id_type vertex_id, const id_type hyperedge_id) noexcept
+    gl_attr_force_inline void bind_head(const id_type vertex_id, const id_type hyperedge_id)
     requires std::same_as<DirectionalTag, hgl::bf_directed_t>
     {
         this->_v_list.bind_head(vertex_id, hyperedge_id);
diff --git a/include/hgl/impl/flat_incidence_matrix.hpp b/include/hgl/impl/flat_incidence_matrix.hpp
index 4a305eb0..b73552f5 100644
--- a/include/hgl/impl/flat_incidence_matrix.hpp
+++ b/include/hgl/impl/flat_incidence_matrix.hpp
@@ -348,21 +348,37 @@ class flat_incidence_matrix<hgl::bf_directed_t, ReprTag> final {
 
     // --- binding methods ---
 
-    gl_attr_force_inline void bind_tail(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_tail(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
+
+        if (this->_matrix[to_idx(major_id), to_idx(minor_id)] == bf_incidence::forward) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the head of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         this->_matrix[to_idx(major_id), to_idx(minor_id)] = bf_incidence::backward;
     }
 
-    gl_attr_force_inline void bind_head(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_head(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
+
+        if (this->_matrix[to_idx(major_id), to_idx(minor_id)] == bf_incidence::backward) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the tail of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         this->_matrix[to_idx(major_id), to_idx(minor_id)] = bf_incidence::forward;
     }
 
-    gl_attr_force_inline void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
+    void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
         this->_matrix[to_idx(major_id), to_idx(minor_id)] = bf_incidence::none;
     }
diff --git a/include/hgl/impl/incidence_list.hpp b/include/hgl/impl/incidence_list.hpp
index a2028475..624339f0 100644
--- a/include/hgl/impl/incidence_list.hpp
+++ b/include/hgl/impl/incidence_list.hpp
@@ -362,25 +362,39 @@ class incidence_list<hgl::bf_directed_t, ReprTag> final {
 
     // --- binding methods ---
 
-    gl_attr_force_inline void bind_tail(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_tail(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
         const auto major_idx = to_idx(major_id);
-        this->_remove_no_align(this->_head_storage[major_idx], minor_id);
+
+        if (this->_contains(this->_head_storage[major_idx], minor_id)) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the head of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         this->_unique_insert(this->_tail_storage[major_idx], minor_id);
     }
 
-    gl_attr_force_inline void bind_head(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_head(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
         const auto major_idx = to_idx(major_id);
-        this->_remove_no_align(this->_tail_storage[major_idx], minor_id);
+
+        if (this->_contains(this->_tail_storage[major_idx], minor_id)) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the tail of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         this->_unique_insert(this->_head_storage[major_idx], minor_id);
     }
 
-    gl_attr_force_inline void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
+    void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
         const auto major_idx = to_idx(major_id);
         this->_remove_no_align(this->_tail_storage[major_idx], minor_id);
@@ -763,14 +777,14 @@ class incidence_list<DirectionalTag, ReprTag> final {
             return this->_e_list.are_bound(vertex_id, hyperedge_id);
     }
 
-    gl_attr_force_inline void bind_tail(const id_type vertex_id, const id_type hyperedge_id) noexcept
+    gl_attr_force_inline void bind_tail(const id_type vertex_id, const id_type hyperedge_id)
     requires std::same_as<DirectionalTag, hgl::bf_directed_t>
     {
         this->_v_list.bind_tail(vertex_id, hyperedge_id);
         this->_e_list.bind_tail(vertex_id, hyperedge_id);
     }
 
-    gl_attr_force_inline void bind_head(const id_type vertex_id, const id_type hyperedge_id) noexcept
+    gl_attr_force_inline void bind_head(const id_type vertex_id, const id_type hyperedge_id)
     requires std::same_as<DirectionalTag, hgl::bf_directed_t>
     {
         this->_v_list.bind_head(vertex_id, hyperedge_id);
diff --git a/include/hgl/impl/incidence_matrix.hpp b/include/hgl/impl/incidence_matrix.hpp
index 3f3c96c6..9179e2b5 100644
--- a/include/hgl/impl/incidence_matrix.hpp
+++ b/include/hgl/impl/incidence_matrix.hpp
@@ -365,21 +365,37 @@ class incidence_matrix<hgl::bf_directed_t, ReprTag> final {
 
     // --- binding methods ---
 
-    gl_attr_force_inline void bind_tail(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_tail(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
+
+        if (this->_matrix[to_idx(major_id)][to_idx(minor_id)] == bf_incidence::forward) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the head of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         this->_matrix[to_idx(major_id)][to_idx(minor_id)] = bf_incidence::backward;
     }
 
-    gl_attr_force_inline void bind_head(
-        const id_type vertex_id, const id_type hyperedge_id
-    ) noexcept {
+    void bind_head(const id_type vertex_id, const id_type hyperedge_id) {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
+
+        if (this->_matrix[to_idx(major_id)][to_idx(minor_id)] == bf_incidence::backward) {
+            throw std::logic_error(std::format(
+                "Tail and head sets must be disjoint: vertex {} is already bound to the tail of "
+                "hyperedge {}.",
+                vertex_id,
+                hyperedge_id
+            ));
+        }
+
         this->_matrix[to_idx(major_id)][to_idx(minor_id)] = bf_incidence::forward;
     }
 
-    gl_attr_force_inline void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
+    void unbind(const id_type vertex_id, const id_type hyperedge_id) noexcept {
         const auto [major_id, minor_id] = layout_tag::majmin(vertex_id, hyperedge_id);
         this->_matrix[to_idx(major_id)][to_idx(minor_id)] = bf_incidence::none;
     }
diff --git a/tests/source/hgl/test_conversion.cpp b/tests/source/hgl/test_conversion.cpp
index 1606a920..eda67cf4 100644
--- a/tests/source/hgl/test_conversion.cpp
+++ b/tests/source/hgl/test_conversion.cpp
@@ -2,7 +2,6 @@
 
 #include <gl/directional_tags.hpp>
 #include <gl/graph_traits.hpp>
-
 #include <hgl/conversion.hpp>
 #include <hgl/hypergraph.hpp>
 #include <hgl/hypergraph_traits.hpp>
diff --git a/tests/source/hgl/test_flat_incidence_list.cpp b/tests/source/hgl/test_flat_incidence_list.cpp
index 4f1e29ae..75744558 100644
--- a/tests/source/hgl/test_flat_incidence_list.cpp
+++ b/tests/source/hgl/test_flat_incidence_list.cpp
@@ -927,7 +927,7 @@ TEST_CASE_FIXTURE(
 
 TEST_CASE_FIXTURE(
     test_bf_directed_vertex_major_flat_incidence_list,
-    "binding methods should rebind the elements if they are already bound"
+    "binding methods should throw if they are already bound"
 ) {
     sut_type sut{constants::n_vertices, constants::n_hyperedges};
     constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
@@ -936,23 +936,33 @@ TEST_CASE_FIXTURE(
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
-    // initial bind
+    // initial bind head
     sut.bind_head(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_head(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
 
     // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
     sut.bind_tail(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_tail(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
     // rebind head
-    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
-    CHECK(sut.is_head(vertex_id, hyperedge_id));
-    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 }
 
 TEST_CASE_FIXTURE(
@@ -1060,32 +1070,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (std::size_t i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (std::size_t i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
@@ -1411,7 +1423,7 @@ TEST_CASE_FIXTURE(
 
 TEST_CASE_FIXTURE(
     test_bf_directed_hyperedge_major_flat_incidence_list,
-    "binding methods should rebind the elements if they are already bound"
+    "binding methods should throw if they are already bound"
 ) {
     sut_type sut{constants::n_vertices, constants::n_hyperedges};
     constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
@@ -1420,23 +1432,33 @@ TEST_CASE_FIXTURE(
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
-    // initial bind
+    // initial bind head
     sut.bind_head(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_head(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
 
     // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
     sut.bind_tail(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_tail(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
     // rebind head
-    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
-    CHECK(sut.is_head(vertex_id, hyperedge_id));
-    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 }
 
 TEST_CASE_FIXTURE(
@@ -1544,32 +1566,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (std::size_t i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (std::size_t i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
diff --git a/tests/source/hgl/test_flat_incidence_matrix.cpp b/tests/source/hgl/test_flat_incidence_matrix.cpp
index e4064362..9e01f67b 100644
--- a/tests/source/hgl/test_flat_incidence_matrix.cpp
+++ b/tests/source/hgl/test_flat_incidence_matrix.cpp
@@ -1013,6 +1013,46 @@ TEST_CASE_FIXTURE(
     CHECK(std::ranges::contains(vertices, constants::id1));
 }
 
+TEST_CASE_FIXTURE(
+    test_bf_directed_vertex_major_flat_incidence_matrix,
+    "binding methods should throw if they are already bound"
+) {
+    sut_type sut{constants::n_vertices, constants::n_hyperedges};
+    constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
+
+    REQUIRE_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // initial bind head
+    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
+    sut.bind_tail(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // rebind head
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+}
+
 TEST_CASE_FIXTURE(
     test_bf_directed_vertex_major_flat_incidence_matrix, "unbind should clear the corresponding bit"
 ) {
@@ -1113,32 +1153,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (auto i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (auto i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
@@ -1468,6 +1510,46 @@ TEST_CASE_FIXTURE(
     CHECK(std::ranges::contains(vertices, constants::id1));
 }
 
+TEST_CASE_FIXTURE(
+    test_bf_directed_hyperedge_major_flat_incidence_matrix,
+    "binding methods should throw if they are already bound"
+) {
+    sut_type sut{constants::n_vertices, constants::n_hyperedges};
+    constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
+
+    REQUIRE_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // initial bind head
+    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
+    sut.bind_tail(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // rebind head
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+}
+
 TEST_CASE_FIXTURE(
     test_bf_directed_hyperedge_major_flat_incidence_matrix,
     "unbind should clear the corresponding bit"
@@ -1569,32 +1651,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (auto i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (auto i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
diff --git a/tests/source/hgl/test_hypergraph.cpp b/tests/source/hgl/test_hypergraph.cpp
index 48626347..03c6abe2 100644
--- a/tests/source/hgl/test_hypergraph.cpp
+++ b/tests/source/hgl/test_hypergraph.cpp
@@ -5,7 +5,6 @@
 #include "testing/hgl/types.hpp"
 
 #include <gl/attributes/diagnostics.hpp>
-
 #include <hgl/directional_tags.hpp>
 #include <hgl/hypergraph.hpp>
 #include <hgl/hypergraph_traits.hpp>
diff --git a/tests/source/hgl/test_incidence_list.cpp b/tests/source/hgl/test_incidence_list.cpp
index 4cd92169..4b80479a 100644
--- a/tests/source/hgl/test_incidence_list.cpp
+++ b/tests/source/hgl/test_incidence_list.cpp
@@ -7,6 +7,7 @@
 
 #include <algorithm>
 #include <ranges>
+#include <stdexcept>
 
 namespace hgl_testing {
 
@@ -912,7 +913,7 @@ TEST_CASE_FIXTURE(
 
 TEST_CASE_FIXTURE(
     test_bf_directed_vertex_major_incidence_list,
-    "binding methods should rebind the elements if they are already bound"
+    "binding methods should throw if they are already bound"
 ) {
     sut_type sut{constants::n_vertices, constants::n_hyperedges};
     constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
@@ -921,23 +922,33 @@ TEST_CASE_FIXTURE(
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
-    // initial bind
+    // initial bind head
     sut.bind_head(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_head(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
 
     // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
     sut.bind_tail(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_tail(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
     // rebind head
-    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
-    CHECK(sut.is_head(vertex_id, hyperedge_id));
-    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 }
 
 TEST_CASE_FIXTURE(
@@ -1045,32 +1056,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (auto i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (auto i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
@@ -1394,7 +1407,7 @@ TEST_CASE_FIXTURE(
 
 TEST_CASE_FIXTURE(
     test_bf_directed_hyperedge_major_incidence_list,
-    "binding methods should rebind the elements if they are already bound"
+    "binding methods should throw if they are already bound"
 ) {
     sut_type sut{constants::n_vertices, constants::n_hyperedges};
     constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
@@ -1403,23 +1416,33 @@ TEST_CASE_FIXTURE(
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
     REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
-    // initial bind
+    // initial bind head
     sut.bind_head(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_head(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
 
     // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
     sut.bind_tail(vertex_id, hyperedge_id);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
     CHECK(sut.is_tail(vertex_id, hyperedge_id));
     CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 
     // rebind head
-    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
     CHECK(sut.are_bound(vertex_id, hyperedge_id));
-    CHECK(sut.is_head(vertex_id, hyperedge_id));
-    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
 }
 
 TEST_CASE_FIXTURE(
@@ -1526,32 +1549,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (auto i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (auto i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
diff --git a/tests/source/hgl/test_incidence_matrix.cpp b/tests/source/hgl/test_incidence_matrix.cpp
index f0edcaae..ca6f179d 100644
--- a/tests/source/hgl/test_incidence_matrix.cpp
+++ b/tests/source/hgl/test_incidence_matrix.cpp
@@ -1006,6 +1006,46 @@ TEST_CASE_FIXTURE(
     CHECK(std::ranges::contains(vertices, constants::id1));
 }
 
+TEST_CASE_FIXTURE(
+    test_bf_directed_vertex_major_incidence_matrix,
+    "binding methods should throw if they are already bound"
+) {
+    sut_type sut{constants::n_vertices, constants::n_hyperedges};
+    constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
+
+    REQUIRE_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // initial bind head
+    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
+    sut.bind_tail(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // rebind head
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+}
+
 TEST_CASE_FIXTURE(
     test_bf_directed_vertex_major_incidence_matrix, "unbind should clear the corresponding bit"
 ) {
@@ -1106,32 +1146,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (auto i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (auto i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }
 
@@ -1460,6 +1502,46 @@ TEST_CASE_FIXTURE(
     CHECK(std::ranges::contains(vertices, constants::id1));
 }
 
+TEST_CASE_FIXTURE(
+    test_bf_directed_hyperedge_major_incidence_matrix,
+    "binding methods should throw if they are already bound"
+) {
+    sut_type sut{constants::n_vertices, constants::n_hyperedges};
+    constexpr auto vertex_id = constants::id1, hyperedge_id = constants::id2;
+
+    REQUIRE_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+    REQUIRE_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // initial bind head
+    sut.bind_head(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // rebind tail
+    CHECK_THROWS_AS(sut.bind_tail(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_head(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_tail(vertex_id, hyperedge_id));
+
+    // unbind
+    sut.unbind(vertex_id, hyperedge_id);
+    CHECK_FALSE(sut.are_bound(vertex_id, hyperedge_id));
+
+    // initial bind tail
+    sut.bind_tail(vertex_id, hyperedge_id);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+
+    // rebind head
+    CHECK_THROWS_AS(sut.bind_head(vertex_id, hyperedge_id), std::logic_error);
+    CHECK(sut.are_bound(vertex_id, hyperedge_id));
+    CHECK(sut.is_tail(vertex_id, hyperedge_id));
+    CHECK_FALSE(sut.is_head(vertex_id, hyperedge_id));
+}
+
 TEST_CASE_FIXTURE(
     test_bf_directed_hyperedge_major_incidence_matrix, "unbind should clear the corresponding bit"
 ) {
@@ -1560,32 +1642,34 @@ TEST_CASE_FIXTURE(
         CHECK(std::ranges::all_of(sut.tail_size_map(n_elements), is_zero));
     }
 
-    // diagonal = tail, everything else is head
-    for (auto i = 0u; i < n_elements; i++) {
-        for (auto j = 0u; j <= i; j++) {
-            if (i == j)
-                sut.bind_tail(i, j);
-            else
-                sut.bind_head(i, j);
+    SUBCASE("mixed bind") {
+        // diagonal = tail, everything else is head
+        for (auto i = 0u; i < n_elements; i++) {
+            for (auto j = 0u; j <= i; j++) {
+                if (i == j)
+                    sut.bind_tail(i, j);
+                else
+                    sut.bind_head(i, j);
+            }
         }
-    }
 
-    const auto deg_map = sut.degree_map(n_elements);
-    const auto out_deg_map = sut.out_degree_map(n_elements);
-    const auto in_deg_map = sut.in_degree_map(n_elements);
+        const auto deg_map = sut.degree_map(n_elements);
+        const auto out_deg_map = sut.out_degree_map(n_elements);
+        const auto in_deg_map = sut.in_degree_map(n_elements);
 
-    const auto esize_map = sut.hyperedge_size_map(n_elements);
-    const auto tsize_map = sut.tail_size_map(n_elements);
-    const auto hsize_map = sut.head_size_map(n_elements);
+        const auto esize_map = sut.hyperedge_size_map(n_elements);
+        const auto tsize_map = sut.tail_size_map(n_elements);
+        const auto hsize_map = sut.head_size_map(n_elements);
 
-    for (auto i = 0uz; i < n_elements; i++) {
-        CHECK_EQ(deg_map[i], i + 1uz);
-        CHECK_EQ(out_deg_map[i], 1uz);
-        CHECK_EQ(in_deg_map[i], i);
+        for (auto i = 0uz; i < n_elements; i++) {
+            CHECK_EQ(deg_map[i], i + 1uz);
+            CHECK_EQ(out_deg_map[i], 1uz);
+            CHECK_EQ(in_deg_map[i], i);
 
-        CHECK_EQ(esize_map[i], n_elements - i);
-        CHECK_EQ(tsize_map[i], 1uz);
-        CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+            CHECK_EQ(esize_map[i], n_elements - i);
+            CHECK_EQ(tsize_map[i], 1uz);
+            CHECK_EQ(hsize_map[i], n_elements - i - 1uz);
+        }
     }
 }