diff --git a/include/boost/graph/louvain_clustering.hpp b/include/boost/graph/louvain_clustering.hpp new file mode 100644 index 000000000..3b59dd54b --- /dev/null +++ b/include/boost/graph/louvain_clustering.hpp @@ -0,0 +1,474 @@ +//======================================================================= +// Copyright 2026 Becheler Code Labs for C++ Alliance +// Authors: Arnaud Becheler +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= +// +// +// Revision History: +// + +#ifndef BOOST_GRAPH_LOUVAIN_CLUSTERING_HPP +#define BOOST_GRAPH_LOUVAIN_CLUSTERING_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost +{ +namespace louvain_detail +{ + +/// @brief Result of graph aggregation operation. +template +struct aggregation_result +{ + Graph graph; + PartitionMap partition; + std::map internal_weights; + std::map> vertex_mapping; +}; + +// Aggregate graph by collapsing communities into super-nodes. +// Edges between communities are preserved with accumulated weights. +template +auto aggregate( + const Graph& g, + const CommunityMap& communities, + const WeightMap& weight +){ + using vertex_descriptor = typename graph_traits::vertex_descriptor; + using edge_descriptor = typename graph_traits::edge_descriptor; + using vertex_iterator = typename graph_traits::vertex_iterator; + using edge_iterator = typename graph_traits::edge_iterator; + using community_type = typename property_traits::value_type; + using weight_type = typename property_traits::value_type; + using edge_property_t = property; + using aggregated_graph_t = adjacency_list; + using result_t = aggregation_result, vertex_descriptor, weight_type>; + + aggregated_graph_t new_g; + vector_property_map new_community_map; + + std::set unique_communities; + std::map comm_to_vertex; + std::map> vertex_to_originals; + + // collect unique communities + vertex_iterator vi, vi_end; + for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) { + unique_communities.insert(get(communities, *vi)); + } + + // create super-nodes with each their own community + for (const community_type& comm : unique_communities) { + vertex_descriptor new_v = add_vertex(new_g); + comm_to_vertex[comm] = new_v; + vertex_to_originals[new_v] = std::set(); + put(new_community_map, new_v, new_v); + } + + // records which original vertices belong to which super-node + for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) { + community_type c = get(communities, *vi); + vertex_descriptor new_v = comm_to_vertex[c]; + vertex_to_originals[new_v].insert(*vi); + } + + // Build edges with accumulated weights + std::map, weight_type> temp_edge_weights; + std::map temp_internal_weights; + + edge_iterator edge_it, edge_end; + for (tie(edge_it, edge_end) = edges(g); edge_it != edge_end; ++edge_it) { + vertex_descriptor u = source(*edge_it, g); + vertex_descriptor v = target(*edge_it, g); + + community_type c_u = get(communities, u); + community_type c_v = get(communities, v); + + vertex_descriptor new_u = comm_to_vertex[c_u]; + vertex_descriptor new_v = comm_to_vertex[c_v]; + + weight_type w = get(weight, *edge_it); + + // vertices in same community = self loop on super node + if (new_u == new_v) { + auto edge_key = std::make_pair(new_u, new_u); + temp_edge_weights[edge_key] += w; + temp_internal_weights[new_u] += w; + } else { + auto edge_key = std::make_pair(std::min(new_u, new_v), std::max(new_u, new_v)); + temp_edge_weights[edge_key] += w; + } + } + + // add edges to connect super nodes + for (const auto& kv : temp_edge_weights) { + edge_descriptor e; + bool inserted; + tie(e, inserted) = add_edge(kv.first.first, kv.first.second, kv.second, new_g); + } + + return result_t{std::move(new_g), std::move(new_community_map), std::move(temp_internal_weights), std::move(vertex_to_originals)}; +} + +// Track hierarchy of aggregation levels for unfolding partitions. +template +struct hierarchy_t +{ + // Each level maps super-nodes to their constituent vertices from the previous level. + using level_t = std::map>; + std::vector levels; + + void push_level(const level_t& mapping) { + levels.push_back(mapping); + } + + void push_level(level_t&& mapping) { + levels.push_back(std::move(mapping)); + } + + std::size_t size() const { + return levels.size(); + } + + bool empty() const { + return levels.empty(); + } + + const level_t& operator[](std::size_t i) const { + return levels[i]; + } + + template + auto unfold(const CommunityMap& final_partition) const + { + assert(!empty()); + + std::map original_partition; + + for (const auto& kv : final_partition) { + std::set current_nodes; + current_nodes.insert(kv.first); + + // From coarse to fine + for (int level = size() - 1; level >= 0; --level) { + std::set next_nodes; + + for (VertexDescriptor node : current_nodes) { + auto it = levels[level].find(node); + assert(it != levels[level].end()); + next_nodes.insert(it->second.begin(), it->second.end()); + } + + current_nodes = std::move(next_nodes); + } + + // Assign all original vertices to community + for (VertexDescriptor original_v : current_nodes) { + original_partition[original_v] = kv.second; + } + } + + return original_partition; + } +}; + +// Create vector of all vertices. +template +auto get_vertex_vector(const Graph& g) +{ + using vertex_descriptor = typename graph_traits::vertex_descriptor; + using vertex_iterator = typename graph_traits::vertex_iterator; + + std::vector vertices_vec; + vertices_vec.reserve(num_vertices(g)); + vertex_iterator vi, vi_end; + for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) { + vertices_vec.push_back(*vi); + } + return vertices_vec; +} + +}} // namespace boost::louvain_detail + +namespace boost +{ + +template +typename property_traits::value_type +louvain_local_optimization( + const Graph& g, + CommunityMap& communities, + const WeightMap& w, + typename property_traits::value_type min_improvement_inner = typename property_traits::value_type(0.0), + unsigned int seed = 0 +) +{ + using community_type = typename property_traits::value_type; + using weight_type = typename property_traits::value_type; + using vertex_descriptor = typename graph_traits::vertex_descriptor; + using vertex_iterator = typename graph_traits::vertex_iterator; + using out_edge_iterator = typename graph_traits::out_edge_iterator; + + std::size_t n = num_vertices(g); + + // Use vector_property_map for O(1) community access + vector_property_map k; + vector_property_map in; + vector_property_map tot; + weight_type m; + + // populates k[c], in[c], tot[c] + weight_type Q = QualityFunction::quality(g, communities, w, k, in, tot, m); + weight_type Q_new = Q; + std::size_t num_moves = 0; + std::size_t pass_number = 0; + bool has_rolled_back = false; + + // Randomize vertex order once + std::mt19937 gen(seed); + std::vector vertex_order = louvain_detail::get_vertex_vector(g); + std::shuffle(vertex_order.begin(), vertex_order.end(), gen); + + // Pre-allocate neighbor buffers + vector_property_map neigh_weight; + std::vector neigh_comm; + neigh_comm.reserve(100); + + // Pre-allocate rollback buffers (only used if needed) + std::vector saved_partition; + vector_property_map saved_in; + vector_property_map saved_tot; + + do + { + Q = Q_new; + num_moves = 0; + pass_number++; + + // Swap 2 random vertices per pass to escape local minima + if (pass_number > 1 && vertex_order.size() >= 2) { + std::uniform_int_distribution dist(0, vertex_order.size() - 1); + std::size_t idx1 = dist(gen); + std::size_t idx2 = dist(gen); + std::swap(vertex_order[idx1], vertex_order[idx2]); + } + + // Lazy save: only save state if we've previously needed to rollback + if (has_rolled_back) { + saved_partition.resize(num_vertices(g)); + saved_in = in; + saved_tot = tot; + vertex_iterator vi_save, vi_save_end; + for (boost::tie(vi_save, vi_save_end) = vertices(g); vi_save != vi_save_end; ++vi_save) { + saved_partition[get(vertex_index, g, *vi_save)] = get(communities, *vi_save); + } + } + + for(auto v : vertex_order) + { + community_type c_old = get(communities, v); + weight_type k_v = get(k, v); + weight_type w_selfloop = 0; + + // Single-pass neighbor aggregation using pre-allocated vector + neigh_comm.clear(); + out_edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = out_edges(v, g); ei != ei_end; ++ei) { + vertex_descriptor neighbor = target(*ei, g); + weight_type edge_w = get(w, *ei); + if (neighbor == v) { + w_selfloop += edge_w; + } else { + community_type c_neighbor = get(communities, neighbor); + if (get(neigh_weight, c_neighbor) == weight_type(0)) { + neigh_comm.push_back(c_neighbor); + } + put(neigh_weight, c_neighbor, get(neigh_weight, c_neighbor) + edge_w); + } + } + + // Remove v from community + weight_type k_v_in_old = get(neigh_weight, c_old); + QualityFunction::remove(in, tot, c_old, k_v, k_v_in_old, w_selfloop); + + // Find best community + community_type c_best = c_old; + weight_type best_gain = 0; + + for(community_type c_neighbor : neigh_comm) + { + weight_type k_v_in_neighbor = get(neigh_weight, c_neighbor); + weight_type gain = QualityFunction::gain(tot, m, c_neighbor, k_v_in_neighbor, k_v); + + if (gain > best_gain) { + best_gain = gain; + c_best = c_neighbor; + } + } + + // Insert v into best community + if (c_best != c_old) { + put(communities, v, c_best); + num_moves++; + } + + weight_type k_v_in_best = get(neigh_weight, c_best); + QualityFunction::insert(in, tot, c_best, k_v, k_v_in_best, w_selfloop); + + // Clear neighbor weights for next vertex + for(community_type c : neigh_comm) { + put(neigh_weight, c, weight_type(0)); + } + } + + // Compute quality from incremental in/tot (no graph traversal) + Q_new = QualityFunction::quality_from_incremental(in, tot, m, n); + + // Rollback if quality didn't improve after moving nodes + // Prevent endless oscillations of vertices: algo gets one extra pass before giving up + if (num_moves > 0 && Q_new <= Q) { + if (has_rolled_back) { + vertex_iterator vi_restore, vi_restore_end; + for (boost::tie(vi_restore, vi_restore_end) = vertices(g); vi_restore != vi_restore_end; ++vi_restore) { + put(communities, *vi_restore, saved_partition[get(vertex_index, g, *vi_restore)]); + } + in = saved_in; + tot = saved_tot; + Q_new = Q; + break; + } else { + has_rolled_back = true; + } + } + + } while (num_moves > 0 && (Q_new - Q) > min_improvement_inner); + + return Q_new; +} + +template +typename property_traits::value_type +louvain_clustering( + const Graph& g0, + ComponentMap components, + const WeightMap& w0, + typename property_traits::value_type min_improvement_inner = typename property_traits::value_type(0.0), + typename property_traits::value_type min_improvement_outer = typename property_traits::value_type(0.0), + unsigned int seed = 0 +){ + using vertex_descriptor = typename graph_traits::vertex_descriptor; + using weight_type = typename property_traits::value_type; + using vertex_iterator = typename graph_traits::vertex_iterator; + + // Initialize each vertex to its own community + vertex_iterator vi, vi_end; + for (boost::tie(vi, vi_end) = vertices(g0); vi != vi_end; ++vi) { + put(components, *vi, *vi); + } + + // Run local optimization + weight_type Q = louvain_local_optimization( + g0, components, w0, min_improvement_inner, seed); + + // Build partition vector from current component map + std::vector partition_vec(num_vertices(g0)); + for (boost::tie(vi, vi_end) = vertices(g0); vi != vi_end; ++vi) { + partition_vec[get(vertex_index, g0, *vi)] = get(components, *vi); + } + + louvain_detail::hierarchy_t hierarchy; + auto partition_map_g0 = make_iterator_property_map(partition_vec.begin(), get(vertex_index, g0)); + auto agg_result = louvain_detail::aggregate(g0, partition_map_g0, w0); + + std::size_t prev_n_vertices = num_vertices(g0); + std::size_t iteration = 0; + + // Track best partition across all levels + std::vector best_partition = partition_vec; + weight_type best_Q = Q; + std::size_t best_level = 0; + + while (true) { + iteration++; + // Check convergence: graph didn't get smaller (no communities merged) + std::size_t n_communities = num_vertices(agg_result.graph); + + if (n_communities >= prev_n_vertices || n_communities == 1) { + break; + } + + hierarchy.push_level(std::move(agg_result.vertex_mapping)); + weight_type Q_old = Q; + weight_type Q_agg = louvain_local_optimization( + agg_result.graph, + agg_result.partition, + get(edge_weight, agg_result.graph), + min_improvement_inner, + seed // Same seed across levels + ); + + // Unfold partition to original graph to compute actual Q + std::map agg_partition_map; + vertex_iterator vi_agg, vi_agg_end; + for (boost::tie(vi_agg, vi_agg_end) = vertices(agg_result.graph); vi_agg != vi_agg_end; ++vi_agg) { + agg_partition_map[*vi_agg] = get(agg_result.partition, *vi_agg); + } + + auto unfolded_map = hierarchy.unfold(agg_partition_map); + vertex_iterator vi_orig, vi_orig_end; + for (boost::tie(vi_orig, vi_orig_end) = vertices(g0); vi_orig != vi_orig_end; ++vi_orig) { + partition_vec[get(vertex_index, g0, *vi_orig)] = unfolded_map[*vi_orig]; + } + + // Compute Q on original graph + auto partition_map_check = make_iterator_property_map(partition_vec.begin(), get(vertex_index, g0)); + Q = QualityFunction::quality(g0, partition_map_check, w0); + + // Track best partition + if (Q > best_Q) { + best_Q = Q; + best_partition = partition_vec; + best_level = iteration; + } + + // Stop if quality did not improve + if (Q - Q_old <= min_improvement_outer) { + break; + } + + prev_n_vertices = n_communities; + agg_result = louvain_detail::aggregate(agg_result.graph, agg_result.partition, get(edge_weight, agg_result.graph)); + } + + // Write best partition to output ComponentMap + partition_vec = best_partition; + Q = best_Q; + + for (boost::tie(vi, vi_end) = vertices(g0); vi != vi_end; ++vi) { + put(components, *vi, partition_vec[get(vertex_index, g0, *vi)]); + } + + return Q; +} + +} // namespace boost + +#endif diff --git a/include/boost/graph/louvain_quality_functions.hpp b/include/boost/graph/louvain_quality_functions.hpp new file mode 100644 index 000000000..f31397877 --- /dev/null +++ b/include/boost/graph/louvain_quality_functions.hpp @@ -0,0 +1,304 @@ +//======================================================================= +// Copyright 2026 Becheler Code Labs for C++ Alliance +// Authors: Arnaud Becheler +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= + +#ifndef BOOST_GRAPH_LOUVAIN_QUALITY_FUNCTIONS_HPP +#define BOOST_GRAPH_LOUVAIN_QUALITY_FUNCTIONS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost +{ + +namespace centrality_detail { + + // Detect if vertex_descriptor is integral (vecS) or pointer-like (listS/setS) + template + struct uses_vector_storage : std::is_integral::vertex_descriptor> {}; + + // Detect if type is hashable, but naive, for BGL integral types are hashable + template + struct is_hashable : std::is_integral {}; + + // Vertex property map selector + template ::value> + struct vertex_pmap_selector; + + // vecS specialization uses vector: get(k,v) O(1), put(k, v, val) O(1), space O(V) pre-allocated array + template + struct vertex_pmap_selector { + using type = vector_property_map; + }; + + // listS/setS specialization uses map for safety: get(k,v) O(log V), put(k, v, val) O(log V), space O(V) pre-allocated array + template + struct vertex_pmap_selector { + using type = associative_property_map::vertex_descriptor, ValueType>>; + }; + + // Community storage selector: picks optimal container based on hashability + template ::value> + struct community_storage_selector; + + // Hashable types use unordered_map: get/put O(1) average + template + struct community_storage_selector { + using type = std::unordered_map; + }; + + // Non-hashable types use map: get/put O(log C) + template + struct community_storage_selector { + using type = std::map; + }; +} + +// Modularity: Q = sum_c [ (L_c/m) - (k_c/2m)^2 ] +// L_c = internal edge weight for community c +// k_c = sum of degrees in community c +// m = total edge weight / 2 +struct newman_and_girvan +{ + + /// Compute modularity quality with provided property maps + template + static inline + typename property_traits::value_type + quality( + const Graph& g, + const CommunityMap& communities, + const WeightMap& weights, + VertexDegreeMap k, + CommunityInMap in, + CommunityTotMap tot, + typename property_traits::value_type& m + ) + { + using community_type = typename property_traits::value_type; + using weight_type = typename property_traits::value_type; + using vertex_descriptor = typename graph_traits::vertex_descriptor; + using edge_iterator = typename graph_traits::edge_iterator; + + // Clear all property maps + m = weight_type(0); + + // Collect all communities and initialize maps + std::set communities_set; + typename graph_traits::vertex_iterator vi, vi_end; + for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) { + put(k, *vi, weight_type(0)); + community_type c = get(communities, *vi); + if (communities_set.insert(c).second) { + // First time seeing this community + put(in, c, weight_type(0)); + put(tot, c, weight_type(0)); + } + } + + edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = edges(g); ei != ei_end; ++ei) + { + vertex_descriptor src = source(*ei, g); + vertex_descriptor trg = target(*ei, g); + + weight_type w = get(weights, *ei); + + community_type c_src = get(communities, src); + community_type c_trg = get(communities, trg); + + if (src == trg) { + // Self-loop counts twice (once per endpoint) + put(k, src, get(k, src) + 2 * w); + put(tot, c_src, get(tot, c_src) + 2 * w); + put(in, c_src, get(in, c_src) + 2 * w); + m += 2 * w; + } else { + // Regular edge + put(k, src, get(k, src) + w); + put(k, trg, get(k, trg) + w); + put(tot, c_src, get(tot, c_src) + w); + put(tot, c_trg, get(tot, c_trg) + w); + m += 2 * w; + + if (c_src == c_trg) { + put(in, c_src, get(in, c_src) + 2 * w); + } + } + } + + // m = (sum of all degrees) / 2 + m /= weight_type(2); + + weight_type two_m = weight_type(2) * m; + + // Empty graphs have zero modularity + if (two_m == weight_type(0)){ + return weight_type(0); + } + + // Q formula: sum_c [ (2*L_c - K_c^2/(2m)) ] / (2m) + weight_type Q(0); + for (const auto& c : communities_set) + { + weight_type K_c = get(tot, c); + weight_type two_L_c = get(in, c); + Q += two_L_c - (K_c * K_c) / two_m; + } + Q /= two_m; + return Q; + } + + /// Compute modularity quality (allocates property maps internally) + template + static inline + typename property_traits::value_type + quality(const Graph& g, const CommunityMap& communities, const WeightMap& weights) + { + using community_type = typename property_traits::value_type; + using weight_type = typename property_traits::value_type; + using vertex_descriptor = typename graph_traits::vertex_descriptor; + using k_map_t = typename centrality_detail::vertex_pmap_selector::type; + using community_storage_t = typename centrality_detail::community_storage_selector::type; + + k_map_t k; + community_storage_t in_map; + community_storage_t tot_map; + auto in = make_assoc_property_map(in_map); + auto tot = make_assoc_property_map(tot_map); + weight_type m; + + return quality(g, communities, weights, k, in, tot, m); + } + + // Incremental updates for local optimization + + /** + * Remove vertex from community. + * @param in Property map: community -> internal edge weights + * @param tot Property map: community -> total edge weights + * @param old_comm Community to remove from + * @param k_v vertex total degree + * @param k_v_in_old Sum of edge weights from vertex to vertices in old_comm + * @param w_selfloop Self-loop weight (default 0) + */ + template + static inline void remove( + CommunityInMap in, + CommunityTotMap tot, + CommunityType old_comm, + WeightType k_v, + WeightType k_v_in_old, + WeightType w_selfloop = WeightType(0) + ) + { + put(in, old_comm, get(in, old_comm) - (2 * k_v_in_old + w_selfloop)); + put(tot, old_comm, get(tot, old_comm) - k_v); + } + + /** + * Insert node into community. + * @param in Property map: community -> internal edge weights + * @param tot Property map: community -> total edge weights + * @param new_comm Community to insert into + * @param k_v Node's total degree + * @param k_v_in_new Sum of edge weights from node to vertices in new_comm + * @param w_selfloop Self-loop weight (default 0) + */ + template + static inline void insert( + CommunityInMap in, + CommunityTotMap tot, + CommunityType new_comm, + WeightType k_v, + WeightType k_v_in_new, + WeightType w_selfloop = WeightType(0) + ) + { + put(in, new_comm, get(in, new_comm) + (2 * k_v_in_new + w_selfloop)); + put(tot, new_comm, get(tot, new_comm) + k_v); + } + + /** + * Compute modularity from incrementally maintained in/tot maps. + * Faster than quality() as it doesn't traverse the graph. + * @param in Property map: community -> internal edge weights + * @param tot Property map: community -> total edge weights + * @param m Total edge weight (half sum of all degrees) + * @param num_communities Number of communities to check + * @return Modularity Q + */ + template + static inline WeightType quality_from_incremental( + CommunityInMap in, + CommunityTotMap tot, + WeightType m, + std::size_t num_communities + ) { + if (m == WeightType(0)) { + return WeightType(0); + } + + WeightType Q = 0; + WeightType two_m = WeightType(2) * m; + + for (std::size_t c = 0; c < num_communities; ++c) { + WeightType K_c = get(tot, c); + if (K_c > WeightType(0)) { + WeightType two_L_c = get(in, c); + Q += two_L_c - (K_c * K_c) / two_m; + } + } + + return Q / two_m; + } + + /** + * Compute modularity gain of moving node to target community. + * @param tot Property map: community -> total edge weights + * @param m Total edge weight (half sum of all degrees) + * @param target_comm Community to evaluate + * @param k_v_in_target Sum of edge weights from node to vertices in target_comm + * @param k_v Node's total degree + * @return Modularity gain + */ + template + static inline WeightType gain( + CommunityTotMap tot, + WeightType m, + CommunityType target_comm, + WeightType k_v_in_target, + WeightType k_v + ) { + // Empty graph check + if (m == WeightType(0)) { + return WeightType(0); + } + + // Gain = k_v_in_target - tot[target] * k_v / (2m) + WeightType tot_target = get(tot, target_comm); + return k_v_in_target - (tot_target * k_v) / (WeightType(2) * m); + } + +}; // newman_and_girvan + +} // namespace boost + +#endif // BOOST_GRAPH_LOUVAIN_QUALITY_FUNCTIONS_HPP diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 11cfbdaae..b8ac87e49 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -144,6 +144,8 @@ alias graph_test_regular : [ run eccentricity.cpp ] [ run clustering_coefficient.cpp ] [ run core_numbers_test.cpp ] + [ run louvain_quality_function_test.cpp ] + [ run louvain_clustering_test.cpp ] [ run read_propmap.cpp ] [ run mcgregor_subgraphs_test.cpp /boost/graph//boost_graph ] [ compile grid_graph_cc.cpp ] diff --git a/test/louvain_clustering_test.cpp b/test/louvain_clustering_test.cpp new file mode 100644 index 000000000..98cd4cd89 --- /dev/null +++ b/test/louvain_clustering_test.cpp @@ -0,0 +1,155 @@ +//======================================================================= +// Copyright 2026 +// Author: Becheler Arnaud +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= + +#include +#include +#include +#include +#include + +using Graph = boost::adjacency_list; +using vertex_descriptor = boost::graph_traits::vertex_descriptor; +using edge_descriptor = boost::graph_traits::edge_descriptor; + +bool approx_equal(double a, double b, double epsilon = 1e-6) { + return std::abs(a - b) < epsilon; +} + +// Edge weight aggregation on barbell graph +void test_aggregation() { + Graph g(6); + add_edge(0, 1, g); + add_edge(1, 2, g); + add_edge(0, 2, g); + add_edge(3, 4, g); + add_edge(4, 5, g); + add_edge(3, 5, g); + add_edge(2, 3, g); + + boost::static_property_map weight_map(1.0); + + std::vector partition = {0, 0, 0, 1, 1, 1}; + auto pmap = boost::make_iterator_property_map(partition.begin(), boost::get(boost::vertex_index, g)); + + auto agg = boost::louvain_detail::aggregate(g, pmap, weight_map); + + BOOST_TEST(boost::num_vertices(agg.graph) == 2); + BOOST_TEST(boost::num_edges(agg.graph) == 3); + + // Find the bridge edge (non-self-loop) + auto wmap = get(boost::edge_weight, agg.graph); + double bridge_weight = 0.0; + for (auto e : boost::make_iterator_range(boost::edges(agg.graph))) { + auto src = boost::source(e, agg.graph); + auto trg = boost::target(e, agg.graph); + if (src != trg) { // Not a self-loop + bridge_weight = get(wmap, e); + } + } + + BOOST_TEST(approx_equal(bridge_weight, 1.0)); +} + +// Ring of cliques benchmark (Blondel et al. 2008) +void test_ring_of_cliques() { + const int num_cliques = 30; + const int clique_size = 5; + const int total_nodes = num_cliques * clique_size; + + Graph g(total_nodes); + + for (int c = 0; c < num_cliques; ++c) { + int base = c * clique_size; + for (int i = 0; i < clique_size; ++i) { + for (int j = i + 1; j < clique_size; ++j) { + add_edge(base + i, base + j, g); + } + } + } + + for (int c = 0; c < num_cliques; ++c) { + int next_c = (c + 1) % num_cliques; + add_edge(c * clique_size + (clique_size - 1), + next_c * clique_size, g); + } + + boost::static_property_map weight_map(1.0); + + std::vector clusters(total_nodes); + auto cluster_map = boost::make_iterator_property_map(clusters.begin(), boost::get(boost::vertex_index, g)); + double Q = boost::louvain_clustering(g, cluster_map, weight_map, 1e-7, 0.0, 42); + + std::set unique_communities(clusters.begin(), clusters.end()); + + BOOST_TEST(Q > 0.80); + BOOST_TEST(Q < 0.95); + BOOST_TEST(unique_communities.size() > 1); + BOOST_TEST(unique_communities.size() <= num_cliques); +} + +// Zachary's karate club +void test_karate_club() { + std::vector> karate_edges = { + {0,1}, {0,2}, {0,3}, {0,4}, {0,5}, {0,6}, {0,7}, {0,8}, {0,10}, {0,11}, {0,12}, {0,13}, {0,17}, {0,19}, {0,21}, {0,31}, + {1,2}, {1,3}, {1,7}, {1,13}, {1,17}, {1,19}, {1,21}, {1,30}, + {2,3}, {2,7}, {2,8}, {2,9}, {2,13}, {2,27}, {2,28}, {2,32}, + {3,7}, {3,12}, {3,13}, + {4,6}, {4,10}, + {5,6}, {5,10}, {5,16}, + {6,16}, + {8,30}, {8,32}, {8,33}, + {9,33}, + {13,33}, + {14,32}, {14,33}, + {15,32}, {15,33}, + {18,32}, {18,33}, + {19,33}, + {20,32}, {20,33}, + {22,32}, {22,33}, + {23,25}, {23,27}, {23,29}, {23,32}, {23,33}, + {24,25}, {24,27}, {24,31}, + {25,31}, + {26,29}, {26,33}, + {27,33}, + {28,31}, {28,33}, + {29,32}, {29,33}, + {30,32}, {30,33}, + {31,32}, {31,33}, + {32,33} + }; + + Graph g(34); + for (const auto& edge : karate_edges) { + add_edge(edge.first, edge.second, g); + } + + boost::static_property_map weight_map(1.0); + + std::vector clusters(34); + auto cluster_map = boost::make_iterator_property_map(clusters.begin(), boost::get(boost::vertex_index, g)); + double Q = boost::louvain_clustering(g, cluster_map, weight_map, 1e-6, 0.0, 42); + + std::set unique_communities(clusters.begin(), clusters.end()); + + std::cout << "Karate Club: Q=" << Q << ", communities=" << unique_communities.size() << "\n"; + + // With seed=42, should get deterministic result + BOOST_TEST(Q > 0.39); + BOOST_TEST(Q < 0.43); + BOOST_TEST(unique_communities.size() >= 3); + BOOST_TEST(unique_communities.size() <= 5); + BOOST_TEST(clusters[0] != clusters[33]); +} + +int main() { + test_aggregation(); + test_ring_of_cliques(); + test_karate_club(); + return boost::report_errors(); +} \ No newline at end of file diff --git a/test/louvain_quality_function_test.cpp b/test/louvain_quality_function_test.cpp new file mode 100644 index 000000000..1510032fa --- /dev/null +++ b/test/louvain_quality_function_test.cpp @@ -0,0 +1,184 @@ +//======================================================================= +// Copyright 2026 +// Author: Becheler Arnaud +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= + +#include +#include +#include +#include + +using Graph = boost::adjacency_list; +using vertex_descriptor = boost::graph_traits::vertex_descriptor; +using edge_descriptor = boost::graph_traits::edge_descriptor; + +bool approx_equal(double a, double b, double epsilon = 1e-6) { + return std::abs(a - b) < epsilon; +} + +// Test modularity with different graph partitions +void test_modularity_different_partitions() { + Graph g(4); + add_edge(0, 1, g); + add_edge(1, 2, g); + add_edge(2, 3, g); + add_edge(3, 0, g); + + boost::static_property_map weight_map(1.0); + + // All in same community : Q = 0 + std::vector partition1 = {0, 0, 0, 0}; + auto pmap1 = boost::make_iterator_property_map(partition1.begin(), boost::get(boost::vertex_index, g)); + double Q1 = boost::newman_and_girvan::quality(g, pmap1, weight_map); + BOOST_TEST(approx_equal(Q1, 0.0)); + + // Two communities (0,1) and (2,3), symetric partition : Q = 0 + std::vector partition2 = {0, 0, 1, 1}; + auto pmap2 = boost::make_iterator_property_map(partition2.begin(), boost::get(boost::vertex_index, g)); + double Q2 = boost::newman_and_girvan::quality(g, pmap2, weight_map); + BOOST_TEST(approx_equal(Q2, 0.0)); + + // All separate communities, many inter-community edges: Q negative + std::vector partition3 = {0, 1, 2, 3}; + auto pmap3 = boost::make_iterator_property_map(partition3.begin(), boost::get(boost::vertex_index, g)); + double Q3 = boost::newman_and_girvan::quality(g, pmap3, weight_map); + BOOST_TEST(Q3 < 0.0); +} + +// Test incremental state operations: remove, insert, gain +void test_state_operations() { + Graph g(4); + add_edge(0, 1, g); + add_edge(1, 2, g); + add_edge(2, 3, g); + + boost::static_property_map weight_map(1.0); + + // Initial partition: two communities {0,1} and {2,3} + std::vector partition = {0, 0, 1, 1}; + auto pmap = boost::make_iterator_property_map(partition.begin(), boost::get(boost::vertex_index, g)); + + // Create property maps + boost::vector_property_map k; + std::map in_map, tot_map; + auto in = boost::make_assoc_property_map(in_map); + auto tot = boost::make_assoc_property_map(tot_map); + double m; + + double Q_initial = boost::newman_and_girvan::quality(g, pmap, weight_map, k, in, tot, m); + + BOOST_TEST(approx_equal(Q_initial, 0.166667, 1e-5)); + + // Test remove operation + vertex_descriptor node_to_move = 1; + vertex_descriptor old_comm = 0; + double k_v = get(k, node_to_move); // degree of vertex 1 is 2 + double k_v_in_old = 1.0; // edge to vertex 0 + + boost::newman_and_girvan::remove(in, tot, old_comm, k_v, k_v_in_old); + + // Community 0 originally had vertices {0,1} with tot = 1+2 = 3 (vertex 0 degree 1, vertex 1 degree 2) + // After removing vertex 1 (degree 2), tot should be 3-2 = 1 + BOOST_TEST(approx_equal(get(tot, old_comm), 1.0)); + + vertex_descriptor new_comm = 1; + double k_v_in_new = 1.0; + boost::newman_and_girvan::insert(in, tot, new_comm, k_v, k_v_in_new); + + // Original community 1 was {2,3} with tot = 2+1 = 3 (vertex 2 : degree 2, vertex 3 : degree 1) + // After adding vertex 1 (degree 2), tot should be 3+2 = 5 + BOOST_TEST(approx_equal(get(tot, new_comm), 5.0)); + + // gain = k_v_in_new - (tot[new_comm] * k_v) / (2*m) + // = 1.0 - (5.0 * 2.0) / (2 * 3.0) = 1.0 - 10.0/6.0 ~= -0.667 + double gain = boost::newman_and_girvan::gain(tot, m, new_comm, k_v_in_new, k_v); + BOOST_TEST(approx_equal(gain, -0.666667, 1e-5)); +} + +// Test weighted graph modularity +void test_modularity_weighted() { + Graph g(3); + auto e1 = add_edge(0, 1, g).first; + auto e2 = add_edge(1, 2, g).first; + auto e3 = add_edge(0, 2, g).first; + + std::map weights; + weights[e1] = 1.0; + weights[e2] = 2.0; + weights[e3] = 1.0; + + auto weight_map = boost::make_assoc_property_map(weights); + std::vector partition = {0, 0, 0}; + auto pmap = boost::make_iterator_property_map(partition.begin(), boost::get(boost::vertex_index, g)); + + // Create property maps + boost::vector_property_map k; + std::map in_map, tot_map; + auto in = boost::make_assoc_property_map(in_map); + auto tot = boost::make_assoc_property_map(tot_map); + double m; + + double Q = boost::newman_and_girvan::quality(g, pmap, weight_map, k, in, tot, m); + + // Complete graph in one community : Q=0 + BOOST_TEST(approx_equal(Q, 0.0)); + BOOST_TEST(approx_equal(m, 4.0)); + BOOST_TEST(approx_equal(get(k, 0), 2.0)); + BOOST_TEST(approx_equal(get(k, 1), 3.0)); + BOOST_TEST(approx_equal(get(k, 2), 3.0)); +} + +// Test modularity with self-loop +// 2 vertices, each with self-loop weight 1, connected by edge weight 2, one community +// Manual calculation (self-loops count twice): +// Self-loops count twice in degree: k[0] = 2*1 + 2 = 4, k[1] = 2*1 + 2 = 4 +// m = (k[0] + k[1]) / 2 = 8/2 = 4, so 2m = 8 +// Community 0: tot = 8, in = 2*1 + 2*2 + 2*1 = 8 +// Q = (1/8) * (8 - 64/8) = 0 +void test_modularity_with_selfloop() { + Graph g(2); + auto e00 = add_edge(0, 0, g).first; + auto e01 = add_edge(0, 1, g).first; + auto e11 = add_edge(1, 1, g).first; + + std::map weights; + weights[e00] = 1.0; + weights[e01] = 2.0; + weights[e11] = 1.0; + + auto weight_map = boost::make_assoc_property_map(weights); + + // All in same community + std::vector partition = {0, 0}; + auto pmap = boost::make_iterator_property_map(partition.begin(), boost::get(boost::vertex_index, g)); + + // Create property maps + boost::vector_property_map k; + std::map in_map, tot_map; + auto in = boost::make_assoc_property_map(in_map); + auto tot = boost::make_assoc_property_map(tot_map); + double m; + + double Q = boost::newman_and_girvan::quality(g, pmap, weight_map, k, in, tot, m); + + double expected_Q = 0.0; + + BOOST_TEST(approx_equal(Q, expected_Q)); + BOOST_TEST(approx_equal(m, 4.0)); + BOOST_TEST(approx_equal(get(k, 0), 4.0)); + BOOST_TEST(approx_equal(get(k, 1), 4.0)); + BOOST_TEST(approx_equal(get(tot, 0), 8.0)); + BOOST_TEST(approx_equal(get(in, 0), 8.0)); +} + +int main() { + test_modularity_different_partitions(); + test_state_operations(); + test_modularity_weighted(); + test_modularity_with_selfloop(); + return boost::report_errors(); +}