@@ -71,59 +71,90 @@ py::list _weighted_triangles_and_degree(py::object G, py::object nodes, py::obje
7171 }
7272
7373 size_t compact_size = adj_vectors.size();
74-
75- #pragma omp parallel
76- {
77- std::vector<char > mark (compact_size, 0 );
78-
79- #pragma omp for schedule(dynamic, 64) nowait
80- for (int idx = 0 ; idx < n; ++idx) {
81- node_t i_id = node_ids[idx];
82- size_t i_compact = node_to_compact_idx[i_id];
83- const auto & nbrs = adj_vectors[i_compact];
84- int m = nbrs.size ();
85-
86- if (m < 2 ) {
87- results[idx] = std::make_tuple (i_id, m, 0.0 );
88- continue ;
89- }
90-
91- // 标记邻居
92- for (node_t v : nbrs) {
93- if (node_to_compact_idx.count (v)) {
94- mark[node_to_compact_idx[v]] = 1 ;
95- }
96- }
97-
98- weight_t weighted_triangles = 0 ;
99- for (int a = 0 ; a < m; ++a) {
100- node_t j_id = nbrs[a];
101- if (!node_to_compact_idx.count (j_id)) continue ;
102-
103- weight_t wij = wt (adj, i_id, j_id, weight_key, max_weight);
104- size_t j_compact = node_to_compact_idx[j_id];
105- const auto & j_nbrs = adj_vectors[j_compact];
106-
107- for (node_t k_id : j_nbrs) {
108- if (k_id > j_id && node_to_compact_idx.count (k_id) && mark[node_to_compact_idx[k_id]]) {
109- weight_t wjk = wt (adj, j_id, k_id, weight_key, max_weight);
110- weight_t wki = wt (adj, k_id, i_id, weight_key, max_weight);
111- weighted_triangles += std::cbrt (wij * wjk * wki);
112- }
113- }
74+
75+ // degree-ordering: compute degree and rank for compact indices
76+ std::vector<size_t> deg(compact_size);
77+ for (size_t i = 0; i < compact_size; ++i) deg[i] = adj_vectors[i].size();
78+ std::vector<size_t> idxs(compact_size);
79+ for (size_t i = 0; i < compact_size; ++i) idxs[i] = i;
80+ std::sort(idxs.begin(), idxs.end(), [&](size_t a, size_t b) {
81+ if (deg[a] != deg[b]) return deg[a] < deg[b];
82+ return a < b;
83+ });
84+ std::vector<size_t> rank(compact_size);
85+ for (size_t i = 0; i < compact_size; ++i) rank[idxs[i]] = i;
86+
87+ // create mapping compact_idx -> original node id for fast lookup
88+ std::vector<node_t> compact_to_node(compact_size);
89+ for (const auto &p : node_to_compact_idx) compact_to_node[p.second] = p.first;
90+
91+ // build oriented adjacency (from lower rank -> higher rank)
92+ std::vector<std::vector<size_t>> out_adj(compact_size);
93+ std::vector<std::vector<weight_t>> out_w(compact_size);
94+ for (size_t u_c = 0; u_c < compact_size; ++u_c) {
95+ const auto &nbrs = adj_vectors[u_c];
96+ node_t u_orig = compact_to_node[u_c];
97+ for (node_t v_id : nbrs) {
98+ auto it = node_to_compact_idx.find(v_id);
99+ if (it == node_to_compact_idx.end()) continue;
100+ size_t v_c = it->second;
101+ if (rank[u_c] < rank[v_c]) {
102+ out_adj[u_c].push_back(v_c);
114103 }
115-
116- // 清除标记
117- for (node_t v : nbrs) {
118- if (node_to_compact_idx.count (v)) {
119- mark[node_to_compact_idx[v]] = 0 ;
104+ }
105+ }
106+
107+ // Precompute weight cbrt per oriented edge using adjacency lookups; also sort neighbor lists
108+ for (size_t u_c = 0; u_c < compact_size; ++u_c) {
109+ auto &nbrs = out_adj[u_c];
110+ std::sort(nbrs.begin(), nbrs.end());
111+ auto &wvec = out_w[u_c];
112+ wvec.reserve(nbrs.size());
113+ node_t u_orig = compact_to_node[u_c];
114+ for (size_t v_c : nbrs) {
115+ node_t v_orig = compact_to_node[v_c];
116+ weight_t w = wt(adj, u_orig, v_orig, weight_key, max_weight);
117+ wvec.push_back(std::cbrt(w));
118+ }
119+ }
120+
121+ // triangle accumulation (each triangle counted once); store contributions per compact node
122+ std::vector<weight_t> tri_contrib(compact_size, 0.0);
123+ for (size_t u_c = 0; u_c < compact_size; ++u_c) {
124+ const auto &out_u = out_adj[u_c];
125+ const auto &w_u = out_w[u_c];
126+ for (size_t iu = 0; iu < out_u.size(); ++iu) {
127+ size_t v_c = out_u[iu];
128+ const auto &out_v = out_adj[v_c];
129+ const auto &w_v = out_w[v_c];
130+ size_t p = 0, q = 0;
131+ while (p < out_u.size() && q < out_v.size()) {
132+ if (out_u[p] < out_v[q]) ++p;
133+ else if (out_v[q] < out_u[p]) ++q;
134+ else {
135+ size_t w_c = out_u[p];
136+ weight_t w_uv = w_u[iu];
137+ weight_t w_uw = w_u[p];
138+ weight_t w_vw = w_v[q];
139+ weight_t contrib = w_uv * w_vw * w_uw; // product of cbrts equals cbrt(product)
140+ tri_contrib[u_c] += contrib;
141+ tri_contrib[v_c] += contrib;
142+ tri_contrib[w_c] += contrib;
143+ ++p; ++q;
120144 }
121145 }
122-
123- results[idx] = std::make_tuple (i_id, m, 2 * weighted_triangles);
124146 }
125147 }
126148
149+ // build results for requested nodes
150+ #pragma omp parallel for schedule(static)
151+ for (int idx = 0; idx < n; ++idx) {
152+ node_t id = node_ids[idx];
153+ size_t c = node_to_compact_idx[id];
154+ int m = adj_vectors[c].size();
155+ results[idx] = std::make_tuple(id, m, 2 * tri_contrib[c]);
156+ }
157+
127158 py::list ret;
128159 for (int i = 0; i < n; ++i) {
129160 ret.append(py::make_tuple(G_.id_to_node[py::cast(std::get<0>(results[i]))],
@@ -173,55 +204,64 @@ py::list _triangles_and_degree(py::object G, py::object nodes = py::none()) {
173204 }
174205
175206 size_t compact_size = adj_vectors.size();
176-
177- #pragma omp parallel
178- {
179- std::vector<char > mark (compact_size, 0 );
180-
181- #pragma omp for schedule(dynamic, 64) nowait
182- for (int idx = 0 ; idx < n; ++idx) {
183- node_t v = node_ids[idx];
184- size_t v_compact = node_to_compact_idx[v];
185- const auto & nbrs = adj_vectors[v_compact];
186- int m = nbrs.size ();
187-
188- if (m < 2 ) {
189- results[idx] = std::make_tuple (v, m, 0 );
190- continue ;
191- }
192-
193- // 标记邻居
194- for (node_t w : nbrs) {
195- if (node_to_compact_idx.count (w)) {
196- mark[node_to_compact_idx[w]] = 1 ;
197- }
198- }
199-
200- int64_t ntriangles = 0 ;
201- for (node_t wa : nbrs) {
202- if (!node_to_compact_idx.count (wa)) continue ;
203-
204- size_t wa_compact = node_to_compact_idx[wa];
205- const auto & wa_nbrs = adj_vectors[wa_compact];
206-
207- for (node_t wb : wa_nbrs) {
208- if (wb > wa && node_to_compact_idx.count (wb) && mark[node_to_compact_idx[wb]]) {
209- ++ntriangles;
210- }
211- }
212- }
213-
214- // 清除标记
215- for (node_t w : nbrs) {
216- if (node_to_compact_idx.count (w)) {
217- mark[node_to_compact_idx[w]] = 0 ;
207+
208+ // degree-ordering and oriented adjacency for unweighted triangle counting
209+ std::vector<size_t> deg(compact_size);
210+ for (size_t i = 0; i < compact_size; ++i) deg[i] = adj_vectors[i].size();
211+ std::vector<size_t> idxs(compact_size);
212+ for (size_t i = 0; i < compact_size; ++i) idxs[i] = i;
213+ std::sort(idxs.begin(), idxs.end(), [&](size_t a, size_t b) {
214+ if (deg[a] != deg[b]) return deg[a] < deg[b];
215+ return a < b;
216+ });
217+ std::vector<size_t> rank(compact_size);
218+ for (size_t i = 0; i < compact_size; ++i) rank[idxs[i]] = i;
219+
220+ std::vector<node_t> compact_to_node(compact_size);
221+ for (const auto &p : node_to_compact_idx) compact_to_node[p.second] = p.first;
222+
223+ std::vector<std::vector<size_t>> out_adj(compact_size);
224+ for (size_t u_c = 0; u_c < compact_size; ++u_c) {
225+ node_t u_orig = compact_to_node[u_c];
226+ for (node_t v_id : adj_vectors[u_c]) {
227+ auto it = node_to_compact_idx.find(v_id);
228+ if (it == node_to_compact_idx.end()) continue;
229+ size_t v_c = it->second;
230+ if (rank[u_c] < rank[v_c]) out_adj[u_c].push_back(v_c);
231+ }
232+ std::sort(out_adj[u_c].begin(), out_adj[u_c].end());
233+ }
234+
235+ std::vector<int64_t> tri_count(compact_size, 0);
236+ for (size_t u_c = 0; u_c < compact_size; ++u_c) {
237+ const auto &out_u = out_adj[u_c];
238+ for (size_t i_u = 0; i_u < out_u.size(); ++i_u) {
239+ size_t v_c = out_u[i_u];
240+ const auto &out_v = out_adj[v_c];
241+ size_t p = 0, q = 0;
242+ while (p < out_u.size() && q < out_v.size()) {
243+ if (out_u[p] < out_v[q]) ++p;
244+ else if (out_v[q] < out_u[p]) ++q;
245+ else {
246+ size_t w_c = out_u[p];
247+ ++tri_count[u_c];
248+ ++tri_count[v_c];
249+ ++tri_count[w_c];
250+ ++p; ++q;
218251 }
219252 }
220-
221- results[idx] = std::make_tuple (v, m, ntriangles);
222253 }
223254 }
224255
256+ // build results for requested nodes
257+ #pragma omp parallel for schedule(static)
258+ for (int idx = 0; idx < n; ++idx) {
259+ node_t v = node_ids[idx];
260+ size_t v_c = node_to_compact_idx[v];
261+ int m = adj_vectors[v_c].size();
262+ results[idx] = std::make_tuple(v, m, tri_count[v_c]);
263+ }
264+
225265 py::list ret;
226266 for (int i = 0; i < n; ++i) {
227267 ret.append(py::make_tuple(G_.id_to_node[py::cast(std::get<0>(results[i]))],
0 commit comments