• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

networkit / networkit / 22305088696

23 Feb 2026 11:56AM UTC coverage: 79.493% (+0.005%) from 79.488%
22305088696

push

github

web-flow
Merge pull request #1393 from Schwarf/coverage/fix_gini_coverage

Add unit test for stats.pyx

29538 of 37158 relevant lines covered (79.49%)

2278086.56 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

73.77
/networkit/cpp/centrality/DynTopHarmonicCloseness.cpp
1
/*
2
 * TopCloseness.cpp
3
 *
4
 *  Created on: 28.02.2018
5
 *      Author: nemes, Eugenio Angriman
6
 */
7

8
#include <cmath>
9
#include <limits>
10
#include <omp.h>
11
#include <span>
12

13
#include <networkit/auxiliary/PrioQueue.hpp>
14
#include <networkit/centrality/DynTopHarmonicCloseness.hpp>
15
#include <networkit/components/ConnectedComponents.hpp>
16
#include <networkit/components/StronglyConnectedComponents.hpp>
17
#include <networkit/distance/AffectedNodes.hpp>
18
#include <networkit/graph/BFS.hpp>
19

20
namespace NetworKit {
21

22
DynTopHarmonicCloseness::DynTopHarmonicCloseness(const Graph &G, count k, bool useBFSbound)
4✔
23
    : G(G), k(k), useBFSbound(useBFSbound),
4✔
24
      allScores(G.upperNodeIdBound(), std::numeric_limits<edgeweight>::max()),
4✔
25
      isExact(G.upperNodeIdBound(), false), isValid(G.upperNodeIdBound(), false),
4✔
26
      cutOff(G.upperNodeIdBound(), std::numeric_limits<edgeweight>::max()),
4✔
27
      exactCutOff(G.upperNodeIdBound(), 0), hasComps(false), component(G.upperNodeIdBound()),
4✔
28
      rOld(G.upperNodeIdBound()) {}
8✔
29

30
DynTopHarmonicCloseness::~DynTopHarmonicCloseness() {
8✔
31
    if (hasComps) {
8✔
32
        delete comps;
4✔
33
    }
34
    if (hasWComps) {
8✔
35
        delete wComps;
4✔
36
    }
37
}
8✔
38

39
std::pair<edgeweight, bool> DynTopHarmonicCloseness::BFScut(node v, edgeweight x, count n, count r,
2,036✔
40
                                                            std::vector<uint8_t> &visited,
41
                                                            std::vector<count> &distances,
42
                                                            std::vector<node> &pred,
43
                                                            count &visitedEdges) {
44

45
    double d = 0;
2,036✔
46
    int64_t gamma = 0, nd = 0;
2,036✔
47

48
    distances[v] = 0;
2,036✔
49
    std::queue<node> Q;
2,036✔
50
    std::queue<node> toReset;
2,036✔
51
    Q.push(v);
2,036✔
52
    toReset.push(v);
2,036✔
53
    ++nd;
2,036✔
54

55
    visited[v] = true;
2,036✔
56

57
    /* Clean up the visited vector for the next run.
58
     * We don't need to do this for the distances and pred vector
59
     * because those are only read after visited[u] is true.
60
     */
61
    auto cleanup = [&]() {
2,036✔
62
        while (!toReset.empty()) {
105,883✔
63
            node w = toReset.front();
103,847✔
64
            toReset.pop();
103,847✔
65
            visited[w] = false;
103,847✔
66
        }
67
    };
4,072✔
68

69
    edgeweight c = 0;
2,036✔
70
    edgeweight ctilde = edgeweight(n - 1);
2,036✔
71

72
    do {
73
        node u = Q.front();
64,199✔
74
        Q.pop();
64,199✔
75

76
        if (distances[u] > d) {
64,199✔
77
            d = d + 1;
4,650✔
78

79
            double d2 = static_cast<double>(int64_t(r) - nd);
4,650✔
80
            ctilde =
4,650✔
81
                c + (static_cast<edgeweight>(gamma) / ((d + 2.0) * (d + 1.0))) + (d2 / (d + 2.0));
4,650✔
82

83
            if (ctilde < x) {
4,650✔
84
                exactCutOff[v] = true;
1,730✔
85
                cutOff[v] = d;
1,730✔
86
                cleanup();
1,730✔
87
                return std::make_pair(ctilde, false);
1,730✔
88
            }
89

90
            gamma = 0;
2,920✔
91
        }
92
        bool cont = true;
62,469✔
93
        G.forNeighborsOf(u, [&](node w) {
62,469✔
94
            if (cont) {
290,485✔
95
                visitedEdges++;
289,332✔
96
                if (!visited[w]) {
289,332✔
97
                    distances[w] = distances[u] + 1;
101,811✔
98

99
                    Q.push(w);
101,811✔
100
                    toReset.push(w);
101,811✔
101
                    visited[w] = true;
101,811✔
102

103
                    c += 1.0 / distances[w];
101,811✔
104

105
                    if (!G.isDirected()) {
101,811✔
106
                        gamma += (G.degree(w) - 1); // notice: only because it's undirected
49,835✔
107
                    } else {
108
                        gamma += G.degree(w);
51,976✔
109
                    }
110

111
                    ++nd;
101,811✔
112
                    pred[w] = u;
101,811✔
113
                } else if (distances[w] > 1 && (pred[u] != w)) {
187,521✔
114
                    ctilde = ctilde - 1.0 / (d + 1.0) + 1.0 / (d + 2.0);
112,353✔
115

116
                    if (ctilde < x) {
112,353✔
117
                        cont = false;
229✔
118
                    }
119
                }
120
            }
121
        });
290,485✔
122
        if (ctilde < x) {
62,469✔
123
            exactCutOff[v] = false;
229✔
124
            cutOff[v] = d;
229✔
125
            cleanup();
229✔
126

127
            return std::make_pair(ctilde, false);
229✔
128
        }
129
    } while (!Q.empty());
62,240✔
130

131
    cleanup();
77✔
132
    exactCutOff[v] = false;
77✔
133
    cutOff[v] = std::numeric_limits<edgeweight>::max();
77✔
134

135
    return std::make_pair(c, true);
77✔
136
}
2,036✔
137

138
void DynTopHarmonicCloseness::BFSbound(node source, std::vector<double> &S2, count *visEdges) {
×
139
    count r = 0;
×
140
    n = G.upperNodeIdBound();
×
141
    std::vector<std::vector<node>> levels(n);
×
142
    std::vector<count> nodesPerLev(n, 0);
×
143
    std::vector<count> sumLevs(n, 0);
×
144
    count nLevs = 0;
×
145
    levels[nLevs].clear();
×
146
    double sum_dist = 0;
×
147
    edgeweight level_bound;
148

149
    auto inverseDistance = [&](edgeweight dist) { return 1.0 / dist; };
×
150

151
    Traversal::BFSfrom(G, source, [&](node u, count dist) {
×
152
        sum_dist += dist > 0 ? inverseDistance(dist) : 0;
×
153

154
        ++r;
×
155
        if (dist > nLevs) {
×
156
            sumLevs[nLevs] += nodesPerLev[nLevs];
×
157
            sumLevs[nLevs + 1] = sumLevs[nLevs];
×
158
            ++nLevs;
×
159
            levels[nLevs].clear();
×
160
        }
161
        levels[nLevs].push_back(u);
×
162
        nodesPerLev[nLevs]++;
×
163
    });
×
164
    sumLevs[nLevs] += nodesPerLev[nLevs];
×
165
    if (G.isDirected()) {
×
166
        (*visEdges) += G.numberOfEdges();
×
167
    } else {
168
        (*visEdges) += 2 * G.numberOfEdges();
×
169
    }
170

171
    S2[source] = sum_dist;
×
172

173
    // we compute the bound for the first level
174
    double closeNodes = 0, farNodes = 0;
×
175
    for (count j = 0; j <= nLevs; j++) {
×
176
        if (j <= 2) {
×
177
            closeNodes += nodesPerLev[j];
×
178
        } else {
179
            farNodes += nodesPerLev[j] * inverseDistance(double(std::abs((double)j - 1.)));
×
180
        }
181
    }
182

183
    level_bound = inverseDistance(2.) * closeNodes + farNodes;
×
184

185
    for (count j = 0; j < levels[1].size(); j++) {
×
186
        node w = levels[1][j];
×
187
        // we subtract 2 not to count the node itself
188
        double bound = (level_bound - inverseDistance(2.)
×
189
                        + (inverseDistance(1.) - inverseDistance(2.)) * G.degree(w));
×
190

191
        if (bound < S2[w] && (!G.isDirected() || component[w] == component[source])) {
×
192
            S2[w] = bound;
×
193
        }
194
    }
195

196
    // now we compute it for the other levels
197
    for (count i = 2; i <= nLevs; i++) {
×
198

199
        level_bound = 0;
×
200
        // TODO: OPTIMIZE?
201
        if (!G.isDirected()) {
×
202
            for (count j = 0; j <= nLevs; j++) {
×
203
                level_bound +=
×
204
                    inverseDistance(std::max(2., double(std::abs((double)j - (double)i))))
×
205
                    * nodesPerLev[j];
×
206
            }
207
        } else {
208
            for (count j = 0; j <= nLevs; j++) {
×
209
                level_bound +=
×
210
                    inverseDistance(std::max(2., (double)j - (double)i)) * nodesPerLev[j];
×
211
            }
212
        }
213

214
        for (count j = 0; j < levels[i].size(); ++j) {
×
215
            node w = levels[i][j];
×
216
            double bound = (level_bound - inverseDistance(2.)
×
217
                            + (inverseDistance(1.) - inverseDistance(2.)) * G.degree(w));
×
218

219
            if (bound < S2[w] && (!G.isDirected() || component[w] == component[source])) {
×
220
                S2[w] = bound;
×
221
            }
222
        }
223
    }
224
}
×
225

226
void DynTopHarmonicCloseness::init() {
4✔
227
    n = G.upperNodeIdBound();
4✔
228

229
    assert(n >= k);
4✔
230

231
    topk.clear();
4✔
232
    topk.resize(k);
4✔
233
    topkScores.clear();
4✔
234
    topkScores.resize(k);
4✔
235

236
    nMinCloseness = 0;
4✔
237
    minCloseness = std::numeric_limits<double>::max();
4✔
238
    trail = 0;
4✔
239
}
4✔
240

241
void DynTopHarmonicCloseness::run() {
4✔
242
    init();
4✔
243
    std::vector<bool> toAnalyze(n, true);
4✔
244
    // We compute the number of reachable nodes (or an upper bound)
245
    // only if we use the algorithm for complex networks.
246
    if (!useBFSbound) {
4✔
247
        computeReachableNodes();
4✔
248
    }
249

250
    // Main priority queue with all nodes in order of decreasing degree
251
    Aux::PrioQueue<edgeweight, node> Q1(n);
4✔
252

253
    G.forNodes([&](node v) { Q1.insert(-(n + G.degree(v)), v); });
2,004✔
254

255
    Aux::PrioQueue<edgeweight, node> top(n);
4✔
256

257
    // protects accesses to all shared variables
258
    omp_lock_t lock;
259
    omp_init_lock(&lock);
4✔
260

261
    edgeweight kth = 0;
4✔
262
#pragma omp parallel
4✔
263
    {
264
        std::vector<uint8_t> visited(n, false);
265
        std::vector<count> distances(n);
266
        std::vector<node> pred(n, 0);
267

268
        std::vector<edgeweight> S(n, std::numeric_limits<edgeweight>::max());
269

270
        while (!Q1.empty()) {
271

272
            omp_set_lock(&lock);
273
            if (Q1.empty()) { // The size of Q1 might have changed.
274
                omp_unset_lock(&lock);
275
                break;
276
            }
277
            std::pair<edgeweight, node> extracted = Q1.extractMin();
278

279
            node v = extracted.second;
280
            toAnalyze[v] = false;
281

282
            omp_unset_lock(&lock);
283

284
            // for networks with large diameters: break if the score of the
285
            // current node is smaller than the k-th highest score
286
            if (useBFSbound && allScores[v] < kth && isValid[v]) {
287
                break;
288
            }
289

290
            if (G.degreeOut(v) == 0) {
291
                omp_set_lock(&lock);
292
                allScores[v] = 0;
293
                isExact[v] = 0;
294
                omp_unset_lock(&lock);
295
            } else if (useBFSbound) {
296
                count visEdges;
297

298
                // Perform a complete BFS from v and obtain upper bounds
299
                // for all nodes in the graph
300
                BFSbound(v, S, &visEdges);
301

302
                omp_set_lock(&lock);
303
                allScores[v] = S[v];
304
                isExact[v] = true;
305
                isValid[v] = true;
306
                omp_unset_lock(&lock);
307

308
                // Update the scores of all nodes with the bounds obtained
309
                // by the complete BFS
310
                G.forNodes([&](node u) {
×
311
                    if (allScores[u] > S[u] && toAnalyze[u]) { // This part must be syncrhonized.
×
312
                        omp_set_lock(&lock);
×
313
                        // Have to check again, because the variables  might have changed
314
                        if (allScores[u] > S[u] && toAnalyze[u]) {
×
315
                            allScores[u] = S[u];
×
316
                            isValid[u] = true;
×
317
                            Q1.remove(-u);
×
318
                            Q1.insert(-allScores[u], u);
×
319
                        }
320
                        omp_unset_lock(&lock);
×
321
                    }
322
                });
×
323
            } else {
324
                count edgeCount = 0;
325

326
                // perform a pruned BFS from v in complex networks
327
                std::pair<edgeweight, bool> c =
328
                    BFScut(v, kth, n, r[v], visited, distances, pred, edgeCount);
329

330
                omp_set_lock(&lock);
331
                allScores[v] = c.first;
332
                isExact[v] = c.second;
333
                isValid[v] = true;
334
                omp_unset_lock(&lock);
335
            }
336

337
            omp_set_lock(&lock);
338
            // Insert v into the list with the k most central nodes if
339
            // its score is larger than the k-th largest value
340
            if (isExact[v] && allScores[v] >= kth) {
341
                top.insert(allScores[v], v);
342
                if (top.size() > k) {
343
                    ++trail;
344
                    if (allScores[v] > kth) {
345
                        if (nMinCloseness == trail) {
346
                            while (top.size() > k) {
347
                                top.extractMin();
348
                            }
349
                            trail = 0;
350
                            nMinCloseness = 1;
351
                            if (k > 1) {
352
                                double pqTop = top.peekMin(0).first;
353
                                double pqNext = pqTop;
354
                                top.forElementsWhile([&]() { return pqTop == pqNext; },
24✔
355
                                                     [&](double curKey, node) {
16✔
356
                                                         pqNext = curKey;
16✔
357
                                                         ++nMinCloseness;
16✔
358
                                                     });
16✔
359
                            }
360
                        }
361
                    } else {
362
                        ++nMinCloseness;
363
                    }
364
                } else if (allScores[v] < minCloseness) {
365
                    minCloseness = allScores[v];
366
                    nMinCloseness = 1;
367
                } else if (allScores[v] == minCloseness) {
368
                    ++nMinCloseness;
369
                }
370
            }
371

372
            // Update the k-th largest value for this thread
373
            if (top.size() >= k) {
374
                std::pair<edgeweight, node> elem = top.extractMin();
375
                kth = elem.first;
376
                top.insert(elem.first, elem.second);
377
                if (nMinCloseness == 1) {
378
                    minCloseness = kth;
379
                }
380
            }
381
            omp_unset_lock(&lock);
382
        }
383
    }
384

385
    if (trail > 0) {
4✔
386
        topk.resize(k + trail);
1✔
387
        topkScores.resize(k + trail);
1✔
388
    }
389

390
    // Store the nodes and their closeness centralities
391
    for (int64_t j = top.size() - 1; j >= 0; --j) {
45✔
392
        std::pair<edgeweight, node> elem = top.extractMin();
41✔
393
        topk[j] = elem.second;
41✔
394
        topkScores[j] = elem.first;
41✔
395
    }
396

397
    for (count j = 0; j < k; ++j) {
44✔
398
        top.insert(topkScores[j], topk[j]);
40✔
399
    }
400

401
    for (count i = 0; i < topk.size() - 1; ++i) {
41✔
402
        count toSort = 1;
37✔
403
        while ((i + toSort) < topk.size() && topkScores[i] == topkScores[i + toSort]) {
37✔
404
            ++toSort;
×
405
        }
406

407
        if (toSort > 1) {
37✔
408
            auto begin = topk.begin() + i;
×
409
            std::sort(begin, begin + toSort);
×
410
            i += toSort - 1;
×
411
        }
412
    }
413

414
    hasRun = true;
4✔
415
}
4✔
416

417
void DynTopHarmonicCloseness::update(GraphEvent event) {
12✔
418

419
    if (event.type == GraphEvent::EDGE_ADDITION) {
12✔
420
        addEdge(event);
8✔
421
    }
422

423
    if (event.type == GraphEvent::EDGE_REMOVAL) {
12✔
424
        removeEdge(event);
4✔
425
    }
426
}
12✔
427

428
void DynTopHarmonicCloseness::addEdge(const GraphEvent &event) {
8✔
429

430
    n = G.upperNodeIdBound();
8✔
431

432
    node eventU = event.u;
8✔
433
    node eventV = event.v;
8✔
434

435
    // Compute the affected nodes
436
    AffectedNodes affectedNodes(G, event);
8✔
437
    affectedNodes.run();
8✔
438

439
    std::vector<node> uniqueAffectedNodes = affectedNodes.getNodes();
8✔
440
    std::vector<edgeweight> distancesFromInsertion = affectedNodes.getDistances();
8✔
441
    std::vector<edgeweight> improvementUpperBounds = affectedNodes.getImprovements();
8✔
442

443
    Aux::PrioQueue<edgeweight, node> Q1(n);
8✔
444

445
    for (node w : uniqueAffectedNodes) {
798✔
446

447
        // Add the improvement bounds for large-diameter networks
448
        if (useBFSbound) {
790✔
449
            isValid[w] = true;
×
450
            isExact[w] = false;
×
451
            allScores[w] += improvementUpperBounds[w];
×
452
        } else {
453
            isValid[w] = false;
790✔
454
        }
455

456
        Q1.insert(-allScores[w], w);
790✔
457
    }
458

459
    allScores[eventU] = affectedNodes.closenessU;
8✔
460
    isExact[eventU] = true;
8✔
461
    isValid[eventU] = true;
8✔
462
    exactCutOff[eventU] = false;
8✔
463
    cutOff[eventU] = std::numeric_limits<edgeweight>::max();
8✔
464

465
    if (!G.isDirected()) {
8✔
466
        isValid[eventV] = true;
4✔
467
        isExact[eventV] = true;
4✔
468
        allScores[eventV] = affectedNodes.closenessV;
4✔
469
        exactCutOff[eventV] = false;
4✔
470
        cutOff[eventV] = std::numeric_limits<edgeweight>::max();
4✔
471
    }
472

473
    Aux::PrioQueue<edgeweight, node> top(n);
8✔
474

475
    count validNodes = 0;
8✔
476
    trail = 0;
8✔
477
    nMinCloseness = 1;
8✔
478
    minCloseness = std::numeric_limits<double>::max();
8✔
479
    // insert all top k nodes with valid closeness into the queue
480
    for (node topNode : topk) {
90✔
481
        if (isValid[topNode] && isExact[topNode]) {
82✔
482
            top.insert(allScores[topNode], topNode);
72✔
483
            if (validNodes == 0) {
72✔
484
                minCloseness = allScores[topNode];
8✔
485
            } else if (minCloseness == allScores[topNode]) {
64✔
486
                ++nMinCloseness;
×
487
            } else {
488
                minCloseness = std::min(minCloseness, allScores[topNode]);
64✔
489
                nMinCloseness = 1;
64✔
490
            }
491
            validNodes++;
72✔
492
            if (validNodes > k) {
72✔
493
                ++trail;
×
494
            }
495
        }
496
    }
497

498
    // Recompute the number of reachable nodes for complex networks
499
    if (!useBFSbound) {
8✔
500
        // Store the old numbers
501
        std::copy(r.begin(), r.end(), rOld.begin());
8✔
502
        updateReachableNodesAfterInsertion(eventU, eventV);
8✔
503
    }
504

505
    // protects accesses to shared data structures
506
    omp_lock_t lock;
507
    omp_init_lock(&lock);
8✔
508

509
    // protects the variables collecting statistics
510
    omp_lock_t statsLock;
511
    omp_init_lock(&statsLock);
8✔
512
    // the new minimum upper bound will be larger or equal to the current
513
    // kth-highest closeness
514
    edgeweight kth = topkScores.back();
8✔
515
#pragma omp parallel
8✔
516
    {
517
        std::vector<uint8_t> visited(n, false);
518
        std::vector<count> distances(n);
519
        count visitedEdges = 0;
520
        std::vector<node> pred(n, 0);
521

522
        std::vector<edgeweight> S(n, std::numeric_limits<edgeweight>::max());
523

524
        while (!Q1.empty()) {
525

526
            omp_set_lock(&lock);
527
            if (Q1.empty()) {
528
                omp_unset_lock(&lock);
529
                break;
530
            }
531
            std::pair<edgeweight, node> extracted = Q1.extractMin();
532
            node v = extracted.second;
533

534
            omp_unset_lock(&lock);
535

536
            if (useBFSbound && allScores[v] < kth && isExact[v]) {
537
                break;
538
            }
539

540
            edgeweight distanceFromInsertedEdge = distancesFromInsertion[v];
541
            edgeweight boundaryUpdateScore =
542
                allScores[v] - 1.0 / (cutOff[v] + 2.0) + 1.0 / (cutOff[v] + 1.0);
543

544
            if ((distanceFromInsertedEdge > cutOff[v] && !isExact[v])
545
                && (!useBFSbound && r[v] <= rOld[v])) {
546
                // our estimate is still valid if the distance of the inserted edge is
547
                // larger than the previous cut-off
548
                isValid[v] = true;
549
                omp_set_lock(&statsLock);
550
                omp_unset_lock(&statsLock);
551
            } else if (distanceFromInsertedEdge == cutOff[v] && !isExact[v]
552
                       && boundaryUpdateScore < kth && (!useBFSbound && r[v] <= rOld[v])) {
553
                // both nodes are affected but not on the same level => cheap update for
554
                // the upper bound
555
                allScores[v] = boundaryUpdateScore;
556
                isValid[v] = true;
557
                omp_set_lock(&statsLock);
558
                omp_unset_lock(&statsLock);
559
            } else if (((!useBFSbound && allScores[v] + improvementUpperBounds[v] < kth)
560
                        || (useBFSbound && allScores[v] < kth))) {
561
                // Adding the improvement bound does not yield a higher upper bound than
562
                // the k-th largest value => add it
563
                if (!useBFSbound) {
564
                    allScores[v] = allScores[v] + improvementUpperBounds[v];
565
                    isExact[v] = false;
566
                    isValid[v] = true;
567
                }
568

569
                omp_set_lock(&statsLock);
570
                omp_unset_lock(&statsLock);
571
            } else if (!(isValid[v] && isExact[v])) {
572
                if (useBFSbound) {
573

574
                    if (isValid[v] && allScores[v] < kth) {
575
                        continue;
576
                    }
577

578
                    count visEdges;
579
                    BFSbound(v, S, &visEdges);
580

581
                    omp_set_lock(&lock);
582
                    allScores[v] = S[v];
583
                    isExact[v] = true;
584
                    isValid[v] = true;
585
                    omp_unset_lock(&lock);
586

587
                    omp_set_lock(&statsLock);
588
                    omp_unset_lock(&statsLock);
589

590
                    G.forNodes([&](node u) {
×
591
                        // This part must be synchronized.
592
                        if ((allScores[u] > S[u] && !isExact[u])) {
×
593
                            omp_set_lock(&lock);
×
594
                            // Have to check again, because the variables might have changed
595
                            if ((allScores[u] > S[u] && !isExact[u])) {
×
596

597
                                allScores[u] = S[u];
×
598
                                isValid[u] = true;
×
599
                                isExact[u] = false;
×
600
                                Q1.remove(-u);
×
601
                                Q1.insert(allScores[u], u);
×
602
                            }
603
                            omp_unset_lock(&lock);
×
604
                        }
605
                    });
×
606
                } else {
607

608
                    std::pair<edgeweight, bool> c =
609
                        BFScut(v, kth, n, r[v], visited, distances, pred, visitedEdges);
610

611
                    allScores[v] = c.first;
612
                    isExact[v] = c.second;
613
                    isValid[v] = true;
614

615
                    omp_set_lock(&statsLock);
616
                    omp_unset_lock(&statsLock);
617
                }
618
            }
619

620
            omp_set_lock(&lock);
621
            if (isExact[v] && allScores[v] >= kth) {
622
                count prevSize = top.size();
623
                top.insert(allScores[v], v);
624
                if (top.size() > std::max(k, prevSize)) {
625
                    ++trail;
626
                    if (allScores[v] > kth) {
627
                        if (nMinCloseness == trail) {
628
                            while (top.size() > k) {
629
                                top.extractMin();
630
                            }
631
                            trail = 0;
632
                            nMinCloseness = 1;
633
                            if (k > 1) {
634
                                double pqTop = top.peekMin(0).first;
635
                                double pqNext = pqTop;
636
                                top.forElementsWhile([&]() { return pqTop == pqNext; },
6✔
637
                                                     [&](double curKey, node) {
4✔
638
                                                         pqNext = curKey;
4✔
639
                                                         ++nMinCloseness;
4✔
640
                                                     });
4✔
641
                            }
642
                        }
643
                    } else {
644
                        ++nMinCloseness;
645
                    }
646
                } else if (allScores[v] < minCloseness) {
647
                    minCloseness = allScores[v];
648
                    nMinCloseness = 1;
649
                } else if (allScores[v] == minCloseness) {
650
                    ++nMinCloseness;
651
                }
652
            }
653

654
            // Update the k-th largest value for this thread
655
            if (top.size() >= k) {
656
                std::pair<edgeweight, node> elem = top.extractMin();
657
                kth = elem.first;
658
                top.insert(elem.first, elem.second);
659
                if (nMinCloseness == 1) {
660
                    minCloseness = kth;
661
                }
662
            }
663
            omp_unset_lock(&lock);
664
        }
665
    }
666

667
    if (trail > 0) {
8✔
668
        topk.resize(k + trail);
×
669
        topkScores.resize(k + trail);
×
670
    }
671

672
    // Store the nodes and their closeness centralities
673
    for (int64_t j = top.size() - 1; j >= 0; --j) {
88✔
674
        std::pair<edgeweight, node> elem = top.extractMin();
80✔
675
        topk[j] = elem.second;
80✔
676
        topkScores[j] = elem.first;
80✔
677
    }
678

679
    for (count i = 0; i < topk.size() - 1; ++i) {
82✔
680
        count toSort = 1;
74✔
681
        while ((i + toSort) < topk.size() && topkScores[i] == topkScores[i + toSort]) {
74✔
682
            ++toSort;
×
683
        }
684
        if (toSort > 1) {
74✔
685
            auto begin = topk.begin() + i;
×
686
            std::sort(begin, begin + toSort);
×
687
            i += toSort - 1;
×
688
        }
689
    }
690
}
8✔
691

692
void DynTopHarmonicCloseness::removeEdge(const GraphEvent &event) {
4✔
693

694
    n = G.upperNodeIdBound();
4✔
695

696
    AffectedNodes affectedNodes(G, event);
4✔
697
    affectedNodes.run();
4✔
698

699
    std::vector<node> uniqueAffectedNodes = affectedNodes.getNodes();
4✔
700

701
    for (node w : uniqueAffectedNodes) {
399✔
702
        isExact[w] = false;
395✔
703
    }
704

705
    Aux::PrioQueue<edgeweight, node> top(n);
4✔
706

707
    // we can abort the update early if none of the top k nodes is actually
708
    // affected
709
    bool hasAffectedTopNode = false;
4✔
710
    for (node v : topk) {
45✔
711
        if (!isExact[v]) {
41✔
712
            hasAffectedTopNode = true;
6✔
713
        }
714
    }
715

716
    if (!hasAffectedTopNode) {
4✔
717
        return;
1✔
718
    }
719

720
    edgeweight kth = 0;
3✔
721
    minCloseness = std::numeric_limits<double>::max();
3✔
722
    nMinCloseness = 1;
3✔
723
    trail = 0;
3✔
724

725
    // This yields to seg. fault when on edge removals on street networks.
726
    if (!useBFSbound) {
3✔
727
        // Since we use BFScut() here, we need the number of reachable nodes in any
728
        // case And now? EA
729
        updateReachableNodesAfterDeletion(event);
3✔
730
    }
731
    Aux::PrioQueue<edgeweight, node> Q(n);
3✔
732

733
    // add all nodes to the queue
734
    G.forNodes([&](node v) { Q.insert(-allScores[v], v); });
1,503✔
735

736
    omp_lock_t lock;
737
    omp_init_lock(&lock);
3✔
738

739
#pragma omp parallel
3✔
740
    {
741
        std::vector<uint8_t> visited(n, false);
742
        std::vector<count> distances(n);
743
        count visitedEdges = 0;
744
        std::vector<node> pred(n, 0);
745

746
        std::vector<edgeweight> S(n, std::numeric_limits<edgeweight>::max());
747
        while (!Q.empty()) {
748

749
            omp_set_lock(&lock);
750
            if (Q.empty()) {
751
                omp_unset_lock(&lock);
752
                break;
753
            }
754
            std::pair<edgeweight, node> elem = Q.extractMin();
755
            node v = elem.second;
756

757
            omp_unset_lock(&lock);
758
            if (allScores[v] < kth) {
759
                break;
760
            }
761

762
            if (!isExact[v] || !isValid[v]) {
763
                if (!useBFSbound) {
764
                    std::pair<edgeweight, bool> c =
765
                        BFScut(v, kth, n, r[v], visited, distances, pred, visitedEdges);
766

767
                    allScores[v] = c.first;
768
                    isExact[v] = c.second;
769
                    isValid[v] = true;
770
                } else {
771
                    BFSbound(v, S, &visitedEdges);
772

773
                    allScores[v] = S[v];
774
                    isExact[v] = true;
775
                    isValid[v] = true;
776
                }
777
            }
778
            omp_set_lock(&lock);
779
            if (isExact[v] && allScores[v] >= kth) {
780
                // Insert node into top queue if closeness is larger than kth
781
                count prevSize = top.size();
782
                top.insert(allScores[v], v);
783
                if (top.size() > std::max(k, prevSize)) {
784
                    ++trail;
785
                    if (allScores[v] > kth) {
786
                        if (nMinCloseness == trail) {
787
                            while (top.size() > k) {
788
                                top.extractMin();
789
                            }
790
                            trail = 0;
791
                            nMinCloseness = 1;
792
                            if (k > 1) {
793
                                double pqTop = top.peekMin(0).first;
794
                                double pqNext = pqTop;
795
                                top.forElementsWhile([&]() { return pqTop == pqNext; },
×
796
                                                     [&](double curKey, node) {
×
797
                                                         pqNext = curKey;
×
798
                                                         ++nMinCloseness;
×
799
                                                     });
×
800
                            }
801
                        }
802
                    } else {
803
                        ++nMinCloseness;
804
                    }
805
                } else if (allScores[v] < minCloseness) {
806
                    minCloseness = allScores[v];
807
                    nMinCloseness = 1;
808
                } else if (allScores[v] == minCloseness) {
809
                    ++nMinCloseness;
810
                }
811
            }
812

813
            // Update kth
814
            if (top.size() >= k) {
815
                std::pair<edgeweight, node> elem = top.extractMin();
816
                kth = elem.first;
817
                top.insert(elem.first, elem.second);
818
                if (nMinCloseness == 1) {
819
                    minCloseness = kth;
820
                }
821
            }
822
            omp_unset_lock(&lock);
823
        }
824
    }
825

826
    // TODO This could go to another method
827
    if (trail > 0) {
3✔
828
        topk.resize(k + trail);
×
829
        topkScores.resize(k + trail);
×
830
    }
831

832
    for (int64_t j = top.size() - 1; j >= 0; j--) {
33✔
833
        std::pair<edgeweight, node> elem = top.extractMin();
30✔
834
        topk[j] = elem.second;
30✔
835
        topkScores[j] = elem.first;
30✔
836
    }
837

838
    for (count j = 0; j < k; j++) {
33✔
839
        top.insert(topkScores[j], topk[j]);
30✔
840
    }
841

842
    for (count i = 0; i < topk.size() - 1; ++i) {
31✔
843
        count toSort = 1;
28✔
844
        while ((i + toSort) < topk.size() && topkScores[i] == topkScores[i + toSort]) {
28✔
845
            ++toSort;
×
846
        }
847
        if (toSort > 1) {
28✔
848
            auto begin = topk.begin() + i;
×
849
            std::sort(begin, begin + toSort);
×
850
            i += toSort - 1;
×
851
        }
852
    }
853
}
6✔
854

855
void DynTopHarmonicCloseness::updateBatch(std::span<const GraphEvent> batch) {
4✔
856
    for (const auto &event : batch) {
8✔
857
        update(event);
4✔
858
    }
859
}
4✔
860

861
void DynTopHarmonicCloseness::reset() {
×
862
    std::fill(isValid.begin(), isValid.end(), false);
×
863
    std::fill(isExact.begin(), isExact.end(), false);
×
864
    std::fill(cutOff.begin(), cutOff.end(), std::numeric_limits<edgeweight>::max());
×
865
}
×
866

867
std::vector<std::pair<node, edgeweight>> DynTopHarmonicCloseness::ranking(bool includeTrail) {
16✔
868
    count nTop = includeTrail ? k + trail : k;
16✔
869
    std::vector<std::pair<node, edgeweight>> ranking(nTop);
16✔
870

871
    for (count i = 0; i < nTop; ++i) {
176✔
872
        ranking[i] = std::make_pair(topk[i], topkScores[i]);
160✔
873
    }
874

875
    return ranking;
16✔
876
}
877

878
void DynTopHarmonicCloseness::computeReachableNodes() {
4✔
879
    if (G.isDirected()) {
4✔
880
        computeReachableNodesDirected();
2✔
881
    } else {
882
        computeReachableNodesUndirected();
2✔
883
    }
884
}
4✔
885

886
void DynTopHarmonicCloseness::computeReachableNodesUndirected() {
2✔
887
    if (!hasComps) {
2✔
888
        comps = new DynConnectedComponents(G);
2✔
889
        hasComps = true;
2✔
890
    }
891

892
    r = std::vector<count>(G.upperNodeIdBound());
2✔
893

894
    comps->run();
2✔
895
    std::map<index, count> sizes = comps->getComponentSizes();
2✔
896
    G.forNodes([&](node v) {
2✔
897
        index cv = comps->componentOfNode(v);
1,000✔
898
        component[v] = cv;
1,000✔
899
        r[v] = sizes[cv];
1,000✔
900
    });
1,000✔
901
}
2✔
902

903
void DynTopHarmonicCloseness::computeReachableNodesDirected() {
2✔
904

905
    r = std::vector<count>(G.upperNodeIdBound());
2✔
906
    wComps = new DynWeaklyConnectedComponents(G);
2✔
907
    wComps->run();
2✔
908
    hasWComps = true;
2✔
909
    std::map<index, count> sizes = wComps->getComponentSizes();
2✔
910
    G.forNodes([&](node w) {
2✔
911
        index cw = wComps->componentOfNode(w);
1,000✔
912
        component[w] = cw;
1,000✔
913
        r[w] = sizes[cw];
1,000✔
914
    });
1,000✔
915
}
2✔
916

917
// TODO: merge in a single method
918
void DynTopHarmonicCloseness::updateReachableNodesAfterInsertion(node u, node v) {
8✔
919

920
    GraphEvent e(GraphEvent::EDGE_ADDITION, u, v);
8✔
921
    if (!G.isDirected()) {
8✔
922

923
        comps->update(e);
4✔
924

925
        std::map<index, count> sizes = comps->getComponentSizes();
4✔
926
        DEBUG(sizes);
4✔
927
        G.forNodes([&](node v) {
4✔
928
            index cv = comps->componentOfNode(v);
2,000✔
929
            component[v] = cv;
2,000✔
930
            r[v] = sizes[cv];
2,000✔
931
        });
2,000✔
932
    } else {
4✔
933
        wComps->update(e);
4✔
934
        // TODO use alias with components so we do not have to replicate the code
935
        std::map<index, count> sizes = wComps->getComponentSizes();
4✔
936
        G.forNodes([&](node) {
4✔
937
            index cv = wComps->componentOfNode(v);
2,000✔
938
            component[v] = cv;
2,000✔
939
            r[v] = sizes[cv];
2,000✔
940
        });
2,000✔
941
    }
4✔
942
}
8✔
943

944
void DynTopHarmonicCloseness::updateReachableNodesAfterDeletion(const GraphEvent &event) {
3✔
945

946
    if (!G.isDirected()) {
3✔
947
        comps->update(event);
2✔
948

949
        std::map<index, count> sizes = comps->getComponentSizes();
2✔
950
        G.forNodes([&](node w) {
2✔
951
            index cv = comps->componentOfNode(w);
1,000✔
952
            component[w] = cv;
1,000✔
953
            r[w] = sizes[cv];
1,000✔
954
        });
1,000✔
955
    } else {
2✔
956
        wComps->update(event);
1✔
957
        std::map<index, count> sizes = wComps->getComponentSizes();
1✔
958
        G.forNodes([&](node w) {
1✔
959
            index cv = wComps->componentOfNode(w);
500✔
960
            component[w] = cv;
500✔
961
            r[w] = sizes[cv];
500✔
962
        });
500✔
963
    }
1✔
964
}
3✔
965

966
} /* namespace NetworKit */
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc