__init__.py 97.2 KB
Newer Older
1
#! /usr/bin/env python
2
# -*- coding: utf-8 -*-
3
#
4
5
# graph_tool -- a general graph manipulation python module
#
Tiago Peixoto's avatar
Tiago Peixoto committed
6
# Copyright (C) 2006-2018 Tiago de Paula Peixoto <tiago@skewed.de>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

21
"""
22
23
``graph_tool.topology`` - Assessing graph topology
--------------------------------------------------
24
25
26
27
28
29
30

Summary
+++++++

.. autosummary::
   :nosignatures:

31
   shortest_distance
Tiago Peixoto's avatar
Tiago Peixoto committed
32
   shortest_path
Tiago Peixoto's avatar
Tiago Peixoto committed
33
34
   all_shortest_paths
   all_predecessors
Tiago Peixoto's avatar
Tiago Peixoto committed
35
   all_paths
Tiago Peixoto's avatar
Tiago Peixoto committed
36
   all_circuits
Tiago Peixoto's avatar
Tiago Peixoto committed
37
   pseudo_diameter
38
   similarity
39
   vertex_similarity
40
   isomorphism
41
42
   subgraph_isomorphism
   mark_subgraph
43
44
   max_cardinality_matching
   max_independent_vertex_set
45
   min_spanning_tree
46
   random_spanning_tree
47
48
49
   dominator_tree
   topological_sort
   transitive_closure
Tiago Peixoto's avatar
Tiago Peixoto committed
50
   tsp_tour
51
   sequential_vertex_coloring
52
53
   label_components
   label_biconnected_components
54
   label_largest_component
55
   label_out_component
56
57
   vertex_percolation
   edge_percolation
Tiago Peixoto's avatar
Tiago Peixoto committed
58
   kcore_decomposition
59
   is_bipartite
Tiago Peixoto's avatar
Tiago Peixoto committed
60
   is_DAG
61
   is_planar
62
   make_maximal_planar
Tiago Peixoto's avatar
Tiago Peixoto committed
63
   edge_reciprocity
64
65
66

Contents
++++++++
67

68
69
"""

70
71
from __future__ import division, absolute_import, print_function

Tiago Peixoto's avatar
Tiago Peixoto committed
72
from .. dl_import import dl_import
73
dl_import("from . import libgraph_tool_topology")
74

75
from .. import _prop, Vector_int32_t, _check_prop_writable, \
76
     _check_prop_scalar, _check_prop_vector, Graph, PropertyMap, GraphView,\
77
     libcore, _get_rng, _degree, perfect_prop_hash, _limit_args
78
from .. stats import label_self_loops
79
import random, sys, numpy, collections
80

81
__all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
82
           "max_cardinality_matching", "max_independent_vertex_set",
83
           "min_spanning_tree", "random_spanning_tree", "dominator_tree",
Tiago Peixoto's avatar
Tiago Peixoto committed
84
           "topological_sort", "transitive_closure", "tsp_tour",
85
86
           "sequential_vertex_coloring", "label_components",
           "label_largest_component", "label_biconnected_components",
87
88
89
90
           "label_out_component", "vertex_percolation", "edge_percolation",
           "kcore_decomposition", "shortest_distance", "shortest_path",
           "all_shortest_paths", "all_predecessors", "all_paths",
           "all_circuits", "pseudo_diameter", "is_bipartite", "is_DAG",
91
92
           "is_planar", "make_maximal_planar", "similarity", "vertex_similarity",
           "edge_reciprocity"]
93

94
def similarity(g1, g2, eweight1=None, eweight2=None, label1=None, label2=None,
95
               norm=True, distance=False, asymmetric=False):
96
97
98
99
100
101
102
    r"""Return the adjacency similarity between the two graphs.

    Parameters
    ----------
    g1 : :class:`~graph_tool.Graph`
        First graph to be compared.
    g2 : :class:`~graph_tool.Graph`
Tiago Peixoto's avatar
Tiago Peixoto committed
103
        Second graph to be compared.
104
105
106
107
    eweight1 : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
        Edge weights for the first graph to be used in comparison.
    eweight2 : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
        Edge weights for the second graph to be used in comparison.
108
109
110
111
112
113
114
115
116
    label1 : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
        Vertex labels for the first graph to be used in comparison. If not
        supplied, the vertex indexes are used.
    label2 : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
        Vertex labels for the second graph to be used in comparison. If not
        supplied, the vertex indexes are used.
    norm : bool (optional, default: ``True``)
        If ``True``, the returned value is normalized by the total number of
        edges.
117
118
119
    distance : bool (optional, default: ``False``)
        If ``True``, the complementary value is returned, i.e. the distance
        between the two graphs.
120
121
122
    asymmetric : bool (optional, default: ``False``)
        If ``True``, the asymmetric similarity of ``g1`` to ``g2`` will be
        computed.
123
124
125
126
127
128
129
130

    Returns
    -------
    similarity : float
        Adjacency similarity value.

    Notes
    -----
Tiago Peixoto's avatar
Tiago Peixoto committed
131
132
133
134
    The adjacency similarity is the sum of equal non-zero entries in the
    adjacency matrix, given a vertex ordering determined by the vertex
    labels. In other words, it counts the number of edges which have the same
    source and target labels in both graphs.
135

136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    More specifically, it is defined as:

    .. math::

       S(\boldsymbol A_1, \boldsymbol A_2) = E - d(\boldsymbol A_1, \boldsymbol A_2)

    where

    .. math::

       d(\boldsymbol A_1, \boldsymbol A_2) = \sum_{i<j} |A_{ij}^{(1)} - A_{ij}^{(2)}|

    is the distance between graphs, and :math:`E=\sum_{i<j}|A_{ij}^{(1)}| +
    |A_{ij}^{(2)}|`.  This definition holds for undirected graphs, otherwise the
    sums go over all directed pairs. If weights are provided, the weighted
    adjacency matrix is used.

    If ``norm == True`` the value returned is
    :math:`S(\boldsymbol A_1, \boldsymbol A_2) / E`.

156
157
158
159
160
    If ``asymmetric == True``, the above is changed so that the comparison is
    made only for entries in :math:`\boldsymbol A_1` that are larger than in :math:`\boldsymbol A_2`, i.e.

    .. math::

161
       d(\boldsymbol A_1, \boldsymbol A_2) = \sum_{i<j} (A_{ij}^{(1)} - A_{ij}^{(2)}) H(A_{ij}^{(1)} - A_{ij}^{(2)}),
162
163
164
165

    where :math:`H(x)` is the unit step function, and the total sum is changed
    accordingly to :math:`E=\sum_{i<j}|A_{ij}^{(1)}|`.

166
167
    The algorithm runs with complexity :math:`O(E_1 + V_1 + E_2 + V_2)`.

168
169
    If enabled during compilation, and the vertex labels are integers bounded by
    the sizes of the graphs, this algorithm runs in parallel.
170

171
172
    Examples
    --------
173
174
175
176
177
178
179
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

180
181
182
183
    >>> g = gt.random_graph(100, lambda: (3,3))
    >>> u = g.copy()
    >>> gt.similarity(u, g)
    1.0
Tiago Peixoto's avatar
Tiago Peixoto committed
184
    >>> gt.random_rewire(u)
Tiago Peixoto's avatar
Tiago Peixoto committed
185
    22
186
    >>> gt.similarity(u, g)
187
    0.04666666666666667
188

189
190
191
192
193
194
    """

    if label1 is None:
        label1 = g1.vertex_index
    if label2 is None:
        label2 = g2.vertex_index
195
196
197
198

    _check_prop_scalar(label1, name="label1")
    _check_prop_scalar(label2, name="label2")

199
    if label1.value_type() != label2.value_type():
200
201
202
203
        try:
            label2 = label2.copy(label1.value_type())
        except ValueError:
            label1 = label1.copy(label2.value_type())
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224

    if eweight1 is None and eweight1 is None:
        ew1 = ew2 = libcore.any()
    else:
        if eweight1 is None:
            eweight1 = g1.new_ep(eweight2.value_type(), 1)
        if eweight2 is None:
            eweight2 = g2.new_ep(eweight1.value_type(), 1)

        _check_prop_scalar(eweight1, name="eweight1")
        _check_prop_scalar(eweight2, name="eweight2")

        if eweight1.value_type() != eweight2.value_type():
            try:
                eweight2 = eweight2.copy(eweight1.value_type())
            except ValueError:
                eweight1 = eweight1.copy(eweight2.value_type())

        ew1 = _prop("e", g1, eweight1)
        ew2 = _prop("e", g2, eweight2)

225
226
    if ((label1.is_writable() and label1.fa.max() > g1.num_vertices()) or
        (label2.is_writable() and label2.fa.max() > g2.num_vertices())):
227
228
        s = libgraph_tool_topology.\
               similarity(g1._Graph__graph, g2._Graph__graph,
229
                          ew1, ew2, _prop("v", g1, label1),
230
                          _prop("v", g2, label2), asymmetric)
Tiago Peixoto's avatar
Tiago Peixoto committed
231
232
233
    else:
        s = libgraph_tool_topology.\
               similarity_fast(g1._Graph__graph, g2._Graph__graph,
234
                               ew1, ew2, _prop("v", g1, label1),
235
                               _prop("v", g2, label2), asymmetric)
236
    if not g1.is_directed() or not g2.is_directed():
237
        s //= 2
238
    if eweight1 is None and eweight1 is None:
239
240
241
242
        if asymmetric:
            E = g1.num_edges()
        else:
            E = g1.num_edges() + g2.num_edges()
243
    else:
244
245
246
247
        if asymmetric:
            E = float(abs(eweight1.fa).sum())
        else:
            E = float(abs(eweight1.fa).sum() + abs(eweight2.fa).sum())
248
249
    if not distance:
        s = E - s
250
    if norm:
251
        return s / E
252
    return s
253

254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
@_limit_args({"sim_type": ["dice", "jaccard", "inv-log-weight"]})
def vertex_similarity(g, sim_type="jaccard", vertex_pairs=None, self_loops=True,
                      sim_map=None):
    r"""Return the similarity between pairs of vertices.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        The graph to be used.
    sim_type : ``str`` (optional, default: ``"jaccard"``)
        Type of similarity to use. This must be one of ``"dice"``, ``"jaccard"``
        or ``"inv-log-weight"``.
    vertex_pairs : iterable of pairs of integers (optional, default: ``None``)
        Pairs of vertices to compute the similarity. If omitted, all pairs will
        be considered.
    self_loops : bool (optional, default: ``True``)
        If ``True``, vertices will be considered adjacent to themselves for the
        purpose of the similarity computation.
    sim_map : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
Tiago Peixoto's avatar
Tiago Peixoto committed
273
        If provided, and ``vertex_pairs is None``, the vertex similarities will
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
        be stored in this vector-valued property. Otherwise, a new one will be
        created.

    Returns
    -------
    similarities : :class:`numpy.ndarray` or :class:`~graph_tool.PropertyMap`
        If ``vertex_pairs`` was supplied, this will be a :class:`numpy.ndarray`
        with the corresponding similarities, otherwise it will be a
        vector-valued vertex :class:`~graph_tool.PropertyMap`, with the
        similarities to all other vertices.

    Notes
    -----
    According to ``sim_type``, this function computes the following similarities:

    ``sim_type == "dice"``

       The Sørensen–Dice similarity [sorensen-dice]_ is twice the number of
292
       common neighbors between two vertices divided by the sum of their
293
       degrees.
294
295
296

    ``sim_type == "jaccard"``

297
298
       The Jaccard similarity [jaccard]_ is the number of common neighbors
       between two vertices divided by the size of the set of all neighbors to
299
300
301
302
303
       both vertices.

    ``sim_type == "inv-log-weight"``

       The inverse log weighted similarity [adamic-friends-2003]_ is the sum of
304
       the weights of common neighbors between two vertices, where the weights
305
       are computed as :math:`1/\log(k)`, with :math:`k` being the degree of the
306
       vertex.
307
308


309
    For directed graphs, only out-neighbors are considered in the above
310
    algorthms (for "inv-log-weight", the in-degrees are used to compute the
311
    weights). To use the in-neighbors instead, a :class:`~graph_tool.GraphView`
312
313
314
315
    should be used to reverse the graph, e.g. ``vertex_similarity(GraphView(g,
    reversed=True))``.

    The algorithm runs with complexity :math:`O(\left<k\right>N^2)` if
Tiago Peixoto's avatar
Tiago Peixoto committed
316
    ``vertex_pairs is None``, otherwise with :math:`O(\left<k\right>P)` where
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
    :math:`P` is the length of ``vertex_pairs``.

    If enabled during compilation, this algorithm runs in parallel.

    Examples
    --------
    .. testcode::
       :hide:

       import matplotlib

    >>> g = gt.collection.data["polbooks"]
    >>> s = gt.vertex_similarity(g, "jaccard")
    >>> color = g.new_vp("double")
    >>> color.a = s[0].a
    >>> gt.graph_draw(g, pos=g.vp.pos, vertex_text=g.vertex_index,
    ...               vertex_color=color, vertex_fill_color=color,
    ...               vcmap=matplotlib.cm.inferno,
    ...               output="polbooks-jaccard.pdf")
    <...>

    .. testcode::
       :hide:

       gt.graph_draw(g, pos=g.vp.pos, vertex_text=g.vertex_index,
                     vertex_color=color, vertex_fill_color=color,
                     vcmap=matplotlib.cm.inferno,
                     output="polbooks-jaccard.png")

    .. figure:: polbooks-jaccard.*
347
       :align: center
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396

       Jaccard similarities to vertex ``0`` in a political books network.

    References
    ----------
    .. [sorensen-dice] https://en.wikipedia.org/wiki/S%C3%B8rensen%E2%80%93Dice_coefficient
    .. [jaccard] https://en.wikipedia.org/wiki/Jaccard_index
    .. [adamic-friends-2003] Lada A. Adamic and Eytan Adar, "Friends and neighbors
       on the Web", Social Networks Volume 25, Issue 3, Pages 211–230 (2003)
       :doi:`10.1016/S0378-8733(03)00009-1`
    .. [liben-nowell-link-prediction-2007] David Liben-Nowell and Jon Kleinberg,
       "The link-prediction problem for social networks", Journal of the
       American Society for Information Science and Technology, Volume 58, Issue
       7, pages 1019–1031 (2007), :doi:`10.1002/asi.20591`
    """

    if vertex_pairs is None:
        if sim_map is None:
            s = g.new_vp("vector<double>")
        else:
            s = sim_map
        if sim_type == "dice":
            libgraph_tool_topology.dice_similarity(g._Graph__graph,
                                                   _prop("v", g, s),
                                                   self_loops)
        elif sim_type == "jaccard":
            libgraph_tool_topology.jaccard_similarity(g._Graph__graph,
                                                      _prop("v", g, s),
                                                      self_loops)
        elif sim_type == "inv-log-weight":
            libgraph_tool_topology.inv_log_weight_similarity(g._Graph__graph,
                                                             _prop("v", g, s))
    else:
        vertex_pairs = numpy.asarray(vertex_pairs, dtype="int64")
        s = numpy.zeros(vertex_pairs.shape[0], dtype="double")
        if sim_type == "dice":
            libgraph_tool_topology.dice_similarity_pairs(g._Graph__graph,
                                                         vertex_pairs,
                                                         s, self_loops)
        elif sim_type == "jaccard":
            libgraph_tool_topology.jaccard_similarity_pairs(g._Graph__graph,
                                                            vertex_pairs,
                                                            s, self_loops)
        elif sim_type == "inv-log-weight":
            libgraph_tool_topology.\
                inv_log_weight_similarity_pairs(g._Graph__graph, vertex_pairs,
                                                s)
    return s

Tiago Peixoto's avatar
Tiago Peixoto committed
397

398
def isomorphism(g1, g2, vertex_inv1=None, vertex_inv2=None, isomap=False):
399
400
    r"""Check whether two graphs are isomorphic.

401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
    Parameters
    ----------
    g1 : :class:`~graph_tool.Graph`
        First graph.
    g2 : :class:`~graph_tool.Graph`
        Second graph.
    vertex_inv1 : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
        Vertex invariant of the first graph. Only vertices with with the same
        invariants are considered in the isomorphism.
    vertex_inv2 : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
        Vertex invariant of the second graph. Only vertices with with the same
        invariants are considered in the isomorphism.
    isomap : ``bool`` (optional, default: ``False``)
        If ``True``, a vertex :class:`~graph_tool.PropertyMap` with the
        isomorphism mapping is returned as well.

    Returns
    -------
    is_isomorphism : ``bool``
        ``True`` if both graphs are isomorphic, otherwise ``False``.
    isomap : :class:`~graph_tool.PropertyMap`
         Isomorphism mapping corresponding to a property map belonging to the
         first graph which maps its vertices to their corresponding vertices of
         the second graph.
425
426
427

    Examples
    --------
428
429
430
431
432
433
434
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

435
436
437
438
439
440
441
442
443
    >>> g = gt.random_graph(100, lambda: (3,3))
    >>> g2 = gt.Graph(g)
    >>> gt.isomorphism(g, g2)
    True
    >>> g.add_edge(g.vertex(0), g.vertex(1))
    <...>
    >>> gt.isomorphism(g, g2)
    False

444
    """
445
    imap = g1.new_vertex_property("int32_t")
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
    if vertex_inv1 is None:
        vertex_inv1 = g1.degree_property_map("total").copy("int64_t")
    else:
        vertex_inv1 = vertex_inv1.copy("int64_t")
        d = g1.degree_property_map("total")
        vertex_inv1.fa += (vertex_inv1.fa.max() + 1) * d.a
    if vertex_inv2 is None:
        vertex_inv2 = g2.degree_property_map("total").copy("int64_t")
    else:
        vertex_inv2 = vertex_inv2.copy("int64_t")
        d = g2.degree_property_map("total")
        vertex_inv2.fa += (vertex_inv2.fa.max() + 1) * d.a

    inv_max = max(vertex_inv1.fa.max(),vertex_inv2.fa.max()) + 1

    l1 = label_self_loops(g1, mark_only=True)
    if l1.fa.max() > 0:
        g1 = GraphView(g1, efilt=1 - l1.fa)

    l2 = label_self_loops(g2, mark_only=True)
    if l2.fa.max() > 0:
        g2 = GraphView(g2, efilt=1 - l2.fa)

469
    iso = libgraph_tool_topology.\
470
           check_isomorphism(g1._Graph__graph, g2._Graph__graph,
471
472
473
                             _prop("v", g1, vertex_inv1),
                             _prop("v", g2, vertex_inv2),
                             inv_max,
Tiago Peixoto's avatar
Tiago Peixoto committed
474
                             _prop("v", g1, imap))
475
476
477
478
479
    if isomap:
        return iso, imap
    else:
        return iso

Tiago Peixoto's avatar
Tiago Peixoto committed
480

481
def subgraph_isomorphism(sub, g, max_n=0, vertex_label=None, edge_label=None,
482
                         induced=False, subgraph=True, generator=False):
483
    r"""Obtain all subgraph isomorphisms of `sub` in `g` (or at most `max_n` subgraphs, if `max_n > 0`).
484

485

Tiago Peixoto's avatar
Tiago Peixoto committed
486
487
488
489
490
491
    Parameters
    ----------
    sub : :class:`~graph_tool.Graph`
        Subgraph for which to be searched.
    g : :class:`~graph_tool.Graph`
        Graph in which the search is performed.
492
    max_n : int (optional, default: ``0``)
Tiago Peixoto's avatar
Tiago Peixoto committed
493
494
        Maximum number of matches to find. If `max_n == 0`, all matches are
        found.
495
    vertex_label : pair of :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
496
        If provided, this should be a pair of :class:`~graph_tool.PropertyMap`
497
498
499
500
        objects, belonging to ``sub`` and ``g`` (in this order), which specify
        vertex labels which should match, in addition to the topological
        isomorphism.
    edge_label : pair of :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
501
        If provided, this should be a pair of :class:`~graph_tool.PropertyMap`
502
503
504
505
506
507
508
        objects, belonging to ``sub`` and ``g`` (in this order), which specify
        edge labels which should match, in addition to the topological
        isomorphism.
    induced : bool (optional, default: ``False``)
        If ``True``, only node-induced subgraphs are found.
    subgraph : bool (optional, default: ``True``)
        If ``False``, all non-subgraph isomorphisms between `sub` and `g` are
509
        found.
510
511
512
513
    generator : bool (optional, default: ``False``)
        If ``True``, a generator will be returned, instead of a list. This is
        useful if the number of isomorphisms is too large to store in memory. If
        ``generator == True``, the option ``max_n`` is ignored.
Tiago Peixoto's avatar
Tiago Peixoto committed
514
515
516

    Returns
    -------
517
518
519
520
    vertex_maps : list (or generator) of :class:`~graph_tool.PropertyMap` objects
        List (or generator) containing vertex property map objects which
        indicate different isomorphism mappings. The property maps vertices in
        `sub` to the corresponding vertex index in `g`.
Tiago Peixoto's avatar
Tiago Peixoto committed
521
522
523

    Notes
    -----
524
525
526
527
528
    The implementation is based on the VF2 algorithm, introduced by Cordella et al.
    [cordella-improved-2001]_ [cordella-subgraph-2004]_. The spatial complexity
    is of order :math:`O(V)`, where :math:`V` is the (maximum) number of vertices
    of the two graphs. Time complexity is :math:`O(V^2)` in the best case and
    :math:`O(V!\times V)` in the worst case.
529
530
531

    Examples
    --------
532
    >>> from numpy.random import poisson
533
534
535
    >>> g = gt.complete_graph(30)
    >>> sub = gt.complete_graph(10)
    >>> vm = gt.subgraph_isomorphism(sub, g, max_n=100)
536
    >>> print(len(vm))
537
    100
538
    >>> for i in range(len(vm)):
539
540
    ...   g.set_vertex_filter(None)
    ...   g.set_edge_filter(None)
541
    ...   vmask, emask = gt.mark_subgraph(g, sub, vm[i])
542
543
    ...   g.set_vertex_filter(vmask)
    ...   g.set_edge_filter(emask)
544
    ...   assert gt.isomorphism(g, sub)
545
546
547
548
    >>> g.set_vertex_filter(None)
    >>> g.set_edge_filter(None)
    >>> ewidth = g.copy_property(emask, value_type="double")
    >>> ewidth.a += 0.5
Tiago Peixoto's avatar
Tiago Peixoto committed
549
550
551
    >>> ewidth.a *= 2
    >>> gt.graph_draw(g, vertex_fill_color=vmask, edge_color=emask,
    ...               edge_pen_width=ewidth, output_size=(200, 200),
552
    ...               output="subgraph-iso-embed.pdf")
553
    <...>
Tiago Peixoto's avatar
Tiago Peixoto committed
554
    >>> gt.graph_draw(sub, output_size=(200, 200), output="subgraph-iso.pdf")
555
556
    <...>

Tiago Peixoto's avatar
Tiago Peixoto committed
557
558
559
560
561
562
563
564
    .. testcode::
       :hide:

       gt.graph_draw(g, vertex_fill_color=vmask, edge_color=emask,
                     edge_pen_width=ewidth, output_size=(200, 200),
                     output="subgraph-iso-embed.png")
       gt.graph_draw(sub, output_size=(200, 200), output="subgraph-iso.png")

Tiago Peixoto's avatar
Tiago Peixoto committed
565
566
    .. image:: subgraph-iso.*
    .. image:: subgraph-iso-embed.*
567

568

Tiago Peixoto's avatar
Tiago Peixoto committed
569
    **Left:** Subgraph searched, **Right:** One isomorphic subgraph found in main graph.
570
571
572

    References
    ----------
573
574
575
576
577
578
    .. [cordella-improved-2001] L. P. Cordella, P. Foggia, C. Sansone, and M. Vento,
       "An improved algorithm for matching large graphs.", 3rd IAPR-TC15 Workshop
       on Graph-based Representations in Pattern Recognition, pp. 149-159, Cuen, 2001.
       http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.101.5342
    .. [cordella-subgraph-2004] L. P. Cordella, P. Foggia, C. Sansone, and M. Vento,
       "A (Sub)Graph Isomorphism Algorithm for Matching Large Graphs.",
Tiago Peixoto's avatar
Tiago Peixoto committed
579
       IEEE Trans. Pattern Anal. Mach. Intell., vol. 26, no. 10, pp. 1367-1372, 2004.
580
581
       :doi:`10.1109/TPAMI.2004.75`
    .. [boost-subgraph-iso] http://www.boost.org/libs/graph/doc/vf2_sub_graph_iso.html
582
    .. [subgraph-isormophism-wikipedia] http://en.wikipedia.org/wiki/Subgraph_isomorphism_problem
583
584

    """
585
586
    if sub.num_vertices() == 0:
        raise ValueError("Cannot search for an empty subgraph.")
587
588
589
590
    if vertex_label is None:
        vertex_label = (None, None)
    elif vertex_label[0].value_type() != vertex_label[1].value_type():
        raise ValueError("Both vertex label property maps must be of the same type!")
591
592
    elif vertex_label[0].value_type() != "int64_t":
        vertex_label = perfect_prop_hash(vertex_label, htype="int64_t")
593

594
595
596
597
    if edge_label is None:
        edge_label = (None, None)
    elif edge_label[0].value_type() != edge_label[1].value_type():
        raise ValueError("Both edge label property maps must be of the same type!")
598
599
    elif edge_label[0].value_type() != "int64_t":
        edge_label = perfect_prop_hash(edge_label, htype="int64_t")
600

601
602
603
604
605
606
607
608
609
    vmaps = libgraph_tool_topology.\
            subgraph_isomorphism(sub._Graph__graph, g._Graph__graph,
                                 _prop("v", sub, vertex_label[0]),
                                 _prop("v", g, vertex_label[1]),
                                 _prop("e", sub, edge_label[0]),
                                 _prop("e", g, edge_label[1]),
                                 max_n, induced, not subgraph,
                                 generator)
    if generator:
610
        return (PropertyMap(vmap, sub, "v") for vmap in vmaps)
611
    else:
612
        return [PropertyMap(vmap, sub, "v") for vmap in vmaps]
613

Tiago Peixoto's avatar
Tiago Peixoto committed
614

615
def mark_subgraph(g, sub, vmap, vmask=None, emask=None):
616
617
618
619
620
621
622
623
624
    r"""
    Mark a given subgraph `sub` on the graph `g`.

    The mapping must be provided by the `vmap` and `emap` parameters,
    which map vertices/edges of `sub` to indexes of the corresponding
    vertices/edges in `g`.

    This returns a vertex and an edge property map, with value type 'bool',
    indicating whether or not a vertex/edge in `g` corresponds to the subgraph
625
    `sub`.
626
    """
627
    if vmask is None:
628
        vmask = g.new_vertex_property("bool")
629
    if emask is None:
630
631
632
633
634
635
636
637
        emask = g.new_edge_property("bool")

    vmask.a = False
    emask.a = False

    for v in sub.vertices():
        w = g.vertex(vmap[v])
        vmask[w] = True
638
        us = set([g.vertex(vmap[x]) for x in v.out_neighbors()])
639

640
        for ew in w.out_edges():
641
642
643
            if ew.target() in us:
                emask[ew] = True

644
    return vmask, emask
645

Tiago Peixoto's avatar
Tiago Peixoto committed
646

647
def min_spanning_tree(g, weights=None, root=None, tree_map=None):
648
649
650
651
652
653
654
    """
    Return the minimum spanning tree of a given graph.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
655
    weights : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
656
657
        The edge weights. If provided, the minimum spanning tree will minimize
        the edge weights.
658
    root : :class:`~graph_tool.Vertex` (optional, default: `None`)
659
        Root of the minimum spanning tree. If this is provided, Prim's algorithm
660
        is used. Otherwise, Kruskal's algorithm is used.
661
    tree_map : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
        If provided, the edge tree map will be written in this property map.

    Returns
    -------
    tree_map : :class:`~graph_tool.PropertyMap`
        Edge property map with mark the tree edges: 1 for tree edge, 0
        otherwise.

    Notes
    -----
    The algorithm runs with :math:`O(E\log E)` complexity, or :math:`O(E\log V)`
    if `root` is specified.

    Examples
    --------
677
678
679
680
681
682
683
684
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

    >>> from numpy.random import random
685
686
687
    >>> g, pos = gt.triangulation(random((400, 2)) * 10, type="delaunay")
    >>> weight = g.new_edge_property("double")
    >>> for e in g.edges():
Tiago Peixoto's avatar
Tiago Peixoto committed
688
    ...    weight[e] = linalg.norm(pos[e.target()].a - pos[e.source()].a)
689
    >>> tree = gt.min_spanning_tree(g, weights=weight)
690
    >>> gt.graph_draw(g, pos=pos, output="triang_orig.pdf")
691
    <...>
692
693
    >>> u = gt.GraphView(g, efilt=tree)
    >>> gt.graph_draw(u, pos=pos, output="triang_min_span_tree.pdf")
694
695
    <...>

Tiago Peixoto's avatar
Tiago Peixoto committed
696
697
698
699
    .. testcode::
       :hide:

       gt.graph_draw(g, pos=pos, output="triang_orig.png")
700
       gt.graph_draw(u, pos=pos, output="triang_min_span_tree.png")
701

702
    .. image:: triang_orig.*
Tiago Peixoto's avatar
Tiago Peixoto committed
703
        :width: 400px
704
    .. image:: triang_min_span_tree.*
Tiago Peixoto's avatar
Tiago Peixoto committed
705
        :width: 400px
706
707

    *Left:* Original graph, *Right:* The minimum spanning tree.
708
709
710
711
712

    References
    ----------
    .. [kruskal-shortest-1956] J. B. Kruskal.  "On the shortest spanning subtree
       of a graph and the traveling salesman problem",  In Proceedings of the
Tiago Peixoto's avatar
Tiago Peixoto committed
713
714
       American Mathematical Society, volume 7, pages 48-50, 1956.
       :doi:`10.1090/S0002-9939-1956-0078686-7`
715
716
717
718
719
    .. [prim-shortest-1957] R. Prim.  "Shortest connection networks and some
       generalizations",  Bell System Technical Journal, 36:1389-1401, 1957.
    .. [boost-mst] http://www.boost.org/libs/graph/doc/graph_theory_review.html#sec:minimum-spanning-tree
    .. [mst-wiki] http://en.wikipedia.org/wiki/Minimum_spanning_tree
    """
720
    if tree_map is None:
721
722
723
724
        tree_map = g.new_edge_property("bool")
    if tree_map.value_type() != "bool":
        raise ValueError("edge property 'tree_map' must be of value type bool.")

725
726
727
728
729
730
731
732
733
734
735
    u = GraphView(g, directed=False)
    if root is None:
        libgraph_tool_topology.\
               get_kruskal_spanning_tree(u._Graph__graph,
                                         _prop("e", g, weights),
                                         _prop("e", g, tree_map))
    else:
        libgraph_tool_topology.\
               get_prim_spanning_tree(u._Graph__graph, int(root),
                                      _prop("e", g, weights),
                                      _prop("e", g, tree_map))
736
    return tree_map
737

Tiago Peixoto's avatar
Tiago Peixoto committed
738

739
def random_spanning_tree(g, weights=None, root=None, tree_map=None):
740
    r"""Return a random spanning tree of a given graph, which can be directed or
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
    undirected.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    weights : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
        The edge weights. If provided, the probability of a particular spanning
        tree being selected is the product of its edge weights.
    root : :class:`~graph_tool.Vertex` (optional, default: `None`)
        Root of the spanning tree. If not provided, it will be selected randomly.
    tree_map : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
        If provided, the edge tree map will be written in this property map.

    Returns
    -------
    tree_map : :class:`~graph_tool.PropertyMap`
        Edge property map with mark the tree edges: 1 for tree edge, 0
        otherwise.

    Notes
    -----
763
764

    The running time for this algorithm is :math:`O(\tau)`, with :math:`\tau`
765
766
767
768
    being the mean hitting time of a random walk on the graph. In the worse case,
    we have :math:`\tau \sim O(V^3)`, with :math:`V` being the number of
    vertices in the graph. However, in much more typical cases (e.g. sparse
    random graphs) the running time is simply :math:`O(V)`.
769
770
771

    Examples
    --------
772
773
774
775
776
777
778
779
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

    >>> from numpy.random import random
780
    >>> g, pos = gt.triangulation(random((400, 2)), type="delaunay")
781
782
783
784
    >>> weight = g.new_edge_property("double")
    >>> for e in g.edges():
    ...    weight[e] = linalg.norm(pos[e.target()].a - pos[e.source()].a)
    >>> tree = gt.random_spanning_tree(g, weights=weight)
785
    >>> tree2 = gt.random_spanning_tree(g, weights=weight)
786
787
    >>> gt.graph_draw(g, pos=pos, output="rtriang_orig.pdf")
    <...>
788
789
790
791
792
    >>> u = gt.GraphView(g, efilt=tree)
    >>> gt.graph_draw(u, pos=pos, output="triang_random_span_tree.pdf")
    <...>
    >>> u2 = gt.GraphView(g, efilt=tree2)
    >>> gt.graph_draw(u2, pos=pos, output="triang_random_span_tree2.pdf")
793
794
    <...>

Tiago Peixoto's avatar
Tiago Peixoto committed
795
796
797
798
    .. testcode::
       :hide:

       gt.graph_draw(g, pos=pos, output="rtriang_orig.png")
799
800
       gt.graph_draw(u, pos=pos, output="triang_random_span_tree.png")
       gt.graph_draw(u2, pos=pos, output="triang_random_span_tree2.png")
801
802

    .. image:: rtriang_orig.*
803
        :width: 300px
804
    .. image:: triang_random_span_tree.*
805
806
807
        :width: 300px
    .. image:: triang_random_span_tree2.*
        :width: 300px
808

809
810
    *Left:* Original graph, *Middle:* A random spanning tree, *Right:* Another
    random spanning tree
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831

    References
    ----------

    .. [wilson-generating-1996] David Bruce Wilson, "Generating random spanning
       trees more quickly than the cover time", Proceedings of the twenty-eighth
       annual ACM symposium on Theory of computing, Pages 296-303, ACM New York,
       1996, :doi:`10.1145/237814.237880`
    .. [boost-rst] http://www.boost.org/libs/graph/doc/random_spanning_tree.html
    """
    if tree_map is None:
        tree_map = g.new_edge_property("bool")
    if tree_map.value_type() != "bool":
        raise ValueError("edge property 'tree_map' must be of value type bool.")

    if root is None:
        root = g.vertex(numpy.random.randint(0, g.num_vertices()),
                        use_index=False)

    # we need to restrict ourselves to the in-component of root
    l = label_out_component(GraphView(g, reversed=True), root)
832
833
834
    u = GraphView(g, vfilt=l)
    if u.num_vertices() != g.num_vertices():
        raise ValueError("There must be a path from all vertices to the root vertex: %d" % int(root) )
835
836
837
838

    libgraph_tool_topology.\
        random_spanning_tree(g._Graph__graph, int(root),
                             _prop("e", g, weights),
839
                             _prop("e", g, tree_map), _get_rng())
840
841
842
    return tree_map


Tiago Peixoto's avatar
Tiago Peixoto committed
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
def dominator_tree(g, root, dom_map=None):
    """Return a vertex property map the dominator vertices for each vertex.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    root : :class:`~graph_tool.Vertex`
        The root vertex.
    dom_map : :class:`~graph_tool.PropertyMap` (optional, default: None)
        If provided, the dominator map will be written in this property map.

    Returns
    -------
    dom_map : :class:`~graph_tool.PropertyMap`
        The dominator map. It contains for each vertex, the index of its
        dominator vertex.

    Notes
    -----
    A vertex u dominates a vertex v, if every path of directed graph from the
    entry to v must go through u.

    The algorithm runs with :math:`O((V+E)\log (V+E))` complexity.

    Examples
    --------
870
871
872
873
874
875
876
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

Tiago Peixoto's avatar
Tiago Peixoto committed
877
878
879
    >>> g = gt.random_graph(100, lambda: (2, 2))
    >>> tree = gt.min_spanning_tree(g)
    >>> g.set_edge_filter(tree)
880
    >>> root = [v for v in g.vertices() if v.in_degree() == 0]
Tiago Peixoto's avatar
Tiago Peixoto committed
881
    >>> dom = gt.dominator_tree(g, root[0])
882
    >>> print(dom.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
883
884
885
886
887
    [ 0  0  0  0  0  0  0 74  0  0  0 97  0  0  0  0  0  0  0  0  0  0  0  0
      0  0  0  0  0 97  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
      0  0  0  0  0  0  0  0  0  0  0  0 64 67  0  0 67  0  0 74  0  0  0  0
     23  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
      0  7  0  0]
Tiago Peixoto's avatar
Tiago Peixoto committed
888
889
890

    References
    ----------
891
    .. [dominator-bgl] http://www.boost.org/libs/graph/doc/lengauer_tarjan_dominator.htm
Tiago Peixoto's avatar
Tiago Peixoto committed
892
893

    """
894
    if dom_map is None:
Tiago Peixoto's avatar
Tiago Peixoto committed
895
896
897
        dom_map = g.new_vertex_property("int32_t")
    if dom_map.value_type() != "int32_t":
        raise ValueError("vertex property 'dom_map' must be of value type" +
898
899
                         " int32_t.")
    if not g.is_directed():
Tiago Peixoto's avatar
Tiago Peixoto committed
900
        raise ValueError("dominator tree requires a directed graph.")
901
    libgraph_tool_topology.\
Tiago Peixoto's avatar
Tiago Peixoto committed
902
903
904
               dominator_tree(g._Graph__graph, int(root),
                              _prop("v", g, dom_map))
    return dom_map
905

Tiago Peixoto's avatar
Tiago Peixoto committed
906

907
def topological_sort(g):
Tiago Peixoto's avatar
Tiago Peixoto committed
908
909
910
911
912
913
914
    """
    Return the topological sort of the given graph. It is returned as an array
    of vertex indexes, in the sort order.

    Notes
    -----
    The topological sort algorithm creates a linear ordering of the vertices
915
    such that if edge (u,v) appears in the graph, then u comes before v in the
Tiago Peixoto's avatar
Tiago Peixoto committed
916
917
918
919
920
921
    ordering. The graph must be a directed acyclic graph (DAG).

    The time complexity is :math:`O(V + E)`.

    Examples
    --------
922
923
924
925
926
927
928
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

Tiago Peixoto's avatar
Tiago Peixoto committed
929
930
931
932
    >>> g = gt.random_graph(30, lambda: (3, 3))
    >>> tree = gt.min_spanning_tree(g)
    >>> g.set_edge_filter(tree)
    >>> sort = gt.topological_sort(g)
933
    >>> print(sort)
Tiago Peixoto's avatar
Tiago Peixoto committed
934
935
    [28 26 29 27 23 22 18 17 16 20 21 15 12 11 10 25 14  9  8  7  5  3  2 24
      4  6  1  0 19 13]
Tiago Peixoto's avatar
Tiago Peixoto committed
936
937
938

    References
    ----------
939
    .. [topological-boost] http://www.boost.org/libs/graph/doc/topological_sort.html
Tiago Peixoto's avatar
Tiago Peixoto committed
940
941
942
943
    .. [topological-wiki] http://en.wikipedia.org/wiki/Topological_sorting

    """

944
    topological_order = Vector_int32_t()
Tiago Peixoto's avatar
Tiago Peixoto committed
945
946
947
948
    is_DAG = libgraph_tool_topology.\
        topological_sort(g._Graph__graph, topological_order)
    if not is_DAG:
        raise ValueError("Graph is not a directed acylic graph (DAG).");
949
    return topological_order.a[::-1].copy()
950

Tiago Peixoto's avatar
Tiago Peixoto committed
951

952
def transitive_closure(g):
Tiago Peixoto's avatar
Tiago Peixoto committed
953
954
955
956
957
958
959
960
961
962
963
964
965
    """Return the transitive closure graph of g.

    Notes
    -----
    The transitive closure of a graph G = (V,E) is a graph G* = (V,E*) such that
    E* contains an edge (u,v) if and only if G contains a path (of at least one
    edge) from u to v. The transitive_closure() function transforms the input
    graph g into the transitive closure graph tc.

    The time complexity (worst-case) is :math:`O(VE)`.

    Examples
    --------
966
967
968
969
970
971
972
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

Tiago Peixoto's avatar
Tiago Peixoto committed
973
974
975
976
977
    >>> g = gt.random_graph(30, lambda: (3, 3))
    >>> tc = gt.transitive_closure(g)

    References
    ----------
978
    .. [transitive-boost] http://www.boost.org/libs/graph/doc/transitive_closure.html
Tiago Peixoto's avatar
Tiago Peixoto committed
979
980
981
982
    .. [transitive-wiki] http://en.wikipedia.org/wiki/Transitive_closure

    """

983
984
985
986
987
988
989
    if not g.is_directed():
        raise ValueError("graph must be directed for transitive closure.")
    tg = Graph()
    libgraph_tool_topology.transitive_closure(g._Graph__graph,
                                              tg._Graph__graph)
    return tg

Tiago Peixoto's avatar
Tiago Peixoto committed
990

991
def label_components(g, vprop=None, directed=None, attractors=False):
992
    """
993
    Label the components to which each vertex in the graph belongs. If the
994
995
    graph is directed, it finds the strongly connected components.

996
997
998
    A property map with the component labels is returned, together with an
    histogram of component labels.

999
1000
    Parameters
    ----------
1001
    g : :class:`~graph_tool.Graph`
1002
        Graph to be used.
1003
    vprop : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
1004
1005
        Vertex property to store the component labels. If none is supplied, one
        is created.
1006
    directed : bool (optional, default: ``None``)
1007
1008
        Treat graph as directed or not, independently of its actual
        directionality.
1009
1010
1011
1012
    attractors : bool (optional, default: ``False``)
        If ``True``, and the graph is directed, an additional array with Boolean
        values is returned, specifying if the strongly connected components are
        attractors or not.
1013
1014
1015

    Returns
    -------
1016
    comp : :class:`~graph_tool.PropertyMap`
1017
        Vertex property map with component labels.
1018
1019
    hist : :class:`~numpy.ndarray`
        Histogram of component labels.
1020
1021
1022
1023
    is_attractor : :class:`~numpy.ndarray`
        A Boolean array specifying if the strongly connected components are
        attractors or not. This returned only if ``attractors == True``, and the
        graph is directed.
1024
1025
1026
1027
1028
1029

    Notes
    -----
    The components are arbitrarily labeled from 0 to N-1, where N is the total
    number of components.

1030
    The algorithm runs in :math:`O(V + E)` time.
1031
1032
1033

    Examples
    --------
1034
1035
1036
1037
1038
1039
    .. testcode::
       :hide:

       numpy.random.seed(43)
       gt.seed_rng(43)

1040
1041
    >>> g = gt.random_graph(100, lambda: (poisson(2), poisson(2)))
    >>> comp, hist, is_attractor = gt.label_components(g, attractors=True)
1042
    >>> print(comp.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
1043
1044
1045
1046
1047
    [ 9  9  9  9 10  1  9 11 12  9  9  9  9  9  9 13  9  9  9  0  9  9 16  9
      9  3  9  9  4 17  9  9 18  9  9 19 20  9  9  9 14  5  9  9  6  9  9  9
     21  9  9  9  9  9  9  9  9  9  9  9  9  9  9  2  9  8  9 22 15  9  9  9
      9  9 23 25  9  9 26 27 28 29 30  9  9  9  9  9  9 31  9  9  9  9  9 32
      9  9  7 24]
1048
    >>> print(hist)
Tiago Peixoto's avatar
Tiago Peixoto committed
1049
1050
    [ 1  1  1  1  1  1  1  1  1 68  1  1  1  1  1  1  1  1  1  1  1  1  1  1
      1  1  1  1  1  1  1  1  1]
1051
    >>> print(is_attractor)
Tiago Peixoto's avatar
Tiago Peixoto committed
1052
1053
1054
    [ True  True  True  True  True  True  True  True  True False  True False
     False False False False False False False False False False False False
     False False False False  True False  True False False]
1055
1056
    """

1057
    if vprop is None:
1058
1059
1060
1061
1062
        vprop = g.new_vertex_property("int32_t")

    _check_prop_writable(vprop, name="vprop")
    _check_prop_scalar(vprop, name="vprop")

1063
1064
    if directed is not None:
        g = GraphView(g, directed=directed)
1065

1066
1067
    hist = libgraph_tool_topology.\
               label_components(g._Graph__graph, _prop("v", g, vprop))
1068
1069
1070
1071
1072
1073
1074
1075
1076

    if attractors and g.is_directed() and directed != False:
        is_attractor = numpy.ones(len(hist), dtype="bool")
        libgraph_tool_topology.\
               label_attractors(g._Graph__graph, _prop("v", g, vprop),
                                is_attractor)
        return vprop, hist, is_attractor
    else:
        return vprop, hist
1077
1078
1079
1080


def label_largest_component(g, directed=None):
    """
1081
1082
    Label the largest component in the graph. If the graph is directed, then the
    largest strongly connected component is labelled.
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096

    A property map with a boolean label is returned.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    directed : bool (optional, default:None)
        Treat graph as directed or not, independently of its actual
        directionality.

    Returns
    -------
    comp : :class:`~graph_tool.PropertyMap`
1097
         Boolean vertex property map which labels the largest component.
1098
1099
1100
1101
1102
1103
1104

    Notes
    -----
    The algorithm runs in :math:`O(V + E)` time.

    Examples
    --------
1105
1106
1107
1108
1109
1110
1111
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

1112
1113
    >>> g = gt.random_graph(100, lambda: poisson(1), directed=False)
    >>> l = gt.label_largest_component(g)
1114
    >>> print(l.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
1115
1116
1117
    [0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0
     0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0]
1118
    >>> u = gt.GraphView(g, vfilt=l)   # extract the largest component as a graph
1119
    >>> print(u.num_vertices())
1120
    18
1121
1122
1123
1124
    """

    label = g.new_vertex_property("bool")
    c, h = label_components(g, directed=directed)
1125
    vfilt, inv = g.get_vertex_filter()
1126
    label.fa = c.fa == h.argmax()
1127
    return label
1128

Tiago Peixoto's avatar
Tiago Peixoto committed
1129

1130
def label_out_component(g, root, label=None):
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
    """
    Label the out-component (or simply the component for undirected graphs) of a
    root vertex.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    root : :class:`~graph_tool.Vertex`
        The root vertex.
1141
1142
1143
    label : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
        If provided, this must be an initialized Boolean vertex property map
        where the out-component will be labeled.
1144
1145
1146

    Returns
    -------
1147
    label : :class:`~graph_tool.PropertyMap`
1148
1149
1150
1151
1152
1153
1154
1155
         Boolean vertex property map which labels the out-component.

    Notes
    -----
    The algorithm runs in :math:`O(V + E)` time.

    Examples
    --------
1156
1157
1158
1159
1160
1161
1162
1163
1164
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

    >>> g = gt.random_graph(100, lambda: poisson(2.2), directed=False)
    >>> l = gt.label_out_component(g, g.vertex(2))
1165
    >>> print(l.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
1166
1167
1168
    [1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1
     1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 0 1 0 1 1 1 0 1
     1 1 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0]
1169
1170
1171

    The in-component can be obtained by reversing the graph.

Tiago Peixoto's avatar
Tiago Peixoto committed
1172
    >>> l = gt.label_out_component(gt.GraphView(g, reversed=True, directed=True),
1173
    ...                            g.vertex(1))
1174
    >>> print(l.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
1175
1176
1177
    [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
1178
1179
    """

1180
1181
1182
1183
1184
    if label is None:
        label = g.new_vertex_property("bool")
    elif label.value_type() != "bool":
        raise ValueError("value type of `label` must be `bool`, not %s" %
                         label.value_type())
1185
1186
1187
1188
1189
1190
    libgraph_tool_topology.\
             label_out_component(g._Graph__graph, int(root),
                                 _prop("v", g, label))
    return label


1191
def label_biconnected_components(g, eprop=None, vprop=None):
1192
1193
1194
1195
    """
    Label the edges of biconnected components, and the vertices which are
    articulation points.

1196
1197
1198
1199
    An edge property map with the component labels is returned, together a
    boolean vertex map marking the articulation points, and an histogram of
    component labels.

1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.

    eprop : :class:`~graph_tool.PropertyMap` (optional, default: None)
        Edge property to label the biconnected components.

    vprop : :class:`~graph_tool.PropertyMap` (optional, default: None)
        Vertex property to mark the articulation points. If none is supplied,
        one is created.


    Returns
    -------
    bicomp : :class:`~graph_tool.PropertyMap`
        Edge property map with the biconnected component labels.
    articulation : :class:`~graph_tool.PropertyMap`
        Boolean vertex property map which has value 1 for each vertex which is
        an articulation point, and zero otherwise.
    nc : int
        Number of biconnected components.

    Notes
    -----

    A connected graph is biconnected if the removal of any single vertex (and
    all edges incident on that vertex) can not disconnect the graph. More
    generally, the biconnected components of a graph are the maximal subsets of
    vertices such that the removal of a vertex from a particular component will
    not disconnect the component. Unlike connected components, vertices may
    belong to multiple biconnected components: those vertices that belong to
    more than one biconnected component are called "articulation points" or,
    equivalently, "cut vertices". Articulation points are vertices whose removal
    would increase the number of connected components in the graph. Thus, a
    graph without articulation points is biconnected. Vertices can be present in
    multiple biconnected components, but each edge can only be contained in a
    single biconnected component.

    The algorithm runs in :math:`O(V + E)` time.

    Examples
    --------
1243
1244
1245
1246
1247
1248
1249
    .. testcode::
       :hide:

       import numpy.random
       numpy.random.seed(42)
       gt.seed_rng(42)

Tiago Peixoto's avatar
Tiago Peixoto committed
1250
    >>> g = gt.random_graph(100, lambda: poisson(2), directed=False)
1251
    >>> comp, art, hist = gt.label_biconnected_components(g)
1252
    >>> print(comp.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
1253
1254
1255
1256
1257
    [26 26 26 26 26 26 26 26 19 25 26 26 23 26 26 26 26  6 26 24 18 26 26 13
     26 26 26 26 26 26 26 26 26 26 26 16 29 26 26 26 26 26 26 15 26 26 26 26
     26  0 26 26 12  2 26 26 26 26 26 26 26 26  9  3 26 28 26 26  8 26  4 26
     26 26 14 26 26 26 26 30 11 26 26 26 20 26 26 27 26 33 26 22 17  7  5 32
     21 26  1 10 31]
1258
    >>> print(art.a)
Tiago Peixoto's avatar
Tiago Peixoto committed
1259
1260
1261
    [1 0 1 1 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0
     1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0
     1 0 0 1 0 0 0 1 1 0 0 1 0 0 1 1 0 0 0 1 0 0 1 0 0 0]
1262
    >>> print(hist)
Tiago Peixoto's avatar
Tiago Peixoto committed
1263
1264
    [ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
      1  1 68  1  1  1  1  1  1  1]
1265

1266
    """
1267

1268
    if vprop is None:
1269
        vprop = g.new_vertex_property("bool")
1270
    if eprop is None:
1271
1272
1273
1274
1275
1276
1277
        eprop = g.new_edge_property("int32_t")

    _check_prop_writable(vprop, name="vprop")
    _check_prop_scalar(vprop, name="vprop")
    _check_prop_writable(eprop, name="eprop")
    _check_prop_scalar(eprop, name="eprop")

1278
1279
    g = GraphView(g, directed=False)
    hist = libgraph_tool_topology.\
1280
1281
             label_biconnected_components(g._Graph__graph, _prop("e", g, eprop),
                                          _prop("v", g, vprop))
1282
    return eprop, vprop, hist
1283

1284
1285
1286
def vertex_percolation(g, vertices, second=False):
    """Compute the size of the largest or second-largest component as vertices
    are (virtually)
1287
1288
1289
1290
1291
1292
1293
1294
    removed from the graph.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    vertices : :class:`numpy.ndarray` or iterable of ints
        List of vertices in reversed order of removal.
1295
1296
    second : bool (optional, default: ``False``)
        If ``True``, the size of the second-largest component will be computed.
1297
1298
1299
1300

    Returns
    -------
    size : :class:`numpy.ndarray`