Merge remote-tracking branch 'refs/remotes/trac/u/gh-tabus/exact_algo…

…rithm_for_diameters_of_large_real_directed_graphs' into t/29309/exact_algorithm_for_diameters_of_large_real_directed_graphs

src/sage/graphs/base/static_sparse_graph.pxd

@@ -36,4 +36,5 @@ cdef uint32_t simple_BFS(short_digraph g,
                          uint32_t *distances,
                          uint32_t *predecessors,
                          uint32_t *waiting_list,
+                         int *mask,
                          bitset_t seen)

src/sage/graphs/base/static_sparse_graph.pyx

@@ -422,6 +422,7 @@ cdef uint32_t simple_BFS(short_digraph g,
                          uint32_t *distances,
                          uint32_t *predecessors,
                          uint32_t *waiting_list,
+                         int *mask,
                          bitset_t seen):
     """
     Perform a breadth first search (BFS) using the same method as in
@@ -452,6 +453,12 @@ cdef uint32_t simple_BFS(short_digraph g,
       assumes that this array has already been allocated. However, there is no
       need to initialize it.
 
+    - ``mask`` -- array of size ``n`` to indicate the vertex to be considered
+      during the BFS search from ``source``. These will be the vertex to 
+      which have the same value as the source. This method assumes that this 
+      array has already been allocated. However, it is possible to pass a 
+      ``NULL`` pointer in which case all vertices will be considered.
+
     - ``seen`` -- bitset of size ``n`` that must be initialized before calling
       this method (i.e., bitset_init(seen, n)). However, there is no need to
       clear it.
@@ -492,7 +499,8 @@ cdef uint32_t simple_BFS(short_digraph g,
 
             # If we notice one of these neighbors is not seen yet, we set its
             # parameters and add it to the queue to be explored later.
-            if not bitset_in(seen, u):
+            if not bitset_in(seen, u) and (mask == NULL or mask[source] == mask[u]):
+
                 distances[u] = distances[v] + 1
                 bitset_add(seen, u)
                 waiting_end += 1
@@ -739,7 +747,7 @@ def tarjan_strongly_connected_components(G):
     cdef list int_to_vertex = list(G)
     cdef short_digraph g
     init_short_digraph(g, G, edge_labelled=False, vertex_list=int_to_vertex)
-    cdef int * scc = <int*> mem.malloc(g.n * sizeof(int))
+    cdef int * scc = <int *> mem.malloc(g.n * sizeof(int))
     sig_on()
     cdef int nscc = tarjan_strongly_connected_components_C(g, scc)
     sig_off()
@@ -752,6 +760,112 @@ def tarjan_strongly_connected_components(G):
         output[scc[i]].append(v)
     return output
 
+cdef uint32_t TYY_diameter_C(short_digraph g):
+    from sage.misc.prandom import randint
+    cdef MemoryAllocator mem = MemoryAllocator()
+
+    cdef uint32_t u, v
+    cdef uint32_t D = 0, n = g.n
+    cdef uint32_t idx, aux_min, aux_max, aux_ecc
+
+    cdef short_digraph rev_g
+    init_reverse(rev_g, g)
+
+    cdef int *scc = <int*> mem.malloc(n * sizeof(int))
+    cdef int nscc = tarjan_strongly_connected_components_C(g, scc)
+    cdef uint32_t *ecc = <uint32_t *> mem.malloc(n * sizeof(uint32_t))
+
+    cdef bitset_t seen
+    bitset_init(seen, n)
+    cdef uint32_t *BFS_order = <uint32_t *> mem.malloc(n * sizeof(uint32_t))
+    cdef uint32_t *distance_from = <uint32_t *> mem.malloc(n * sizeof(uint32_t))
+    cdef uint32_t *distance_back = <uint32_t *> mem.malloc(n * sizeof(uint32_t))
+
+    # Initialize distances
+    for idx in range(n):
+        ecc[idx] = UINT32_MAX-1 # so that the first loop will not overflow
+
+    # DoubleSweep
+    v = randint(0, n-1)
+    simple_BFS(g, v, distance_from, NULL, BFS_order, NULL, seen)
+
+    cdef uint32_t scc_sz = 0
+    for u in range(n):
+        if scc[v] == scc[u]:
+            scc_sz += 1
+
+    u = BFS_order[scc_sz-1]
+    simple_BFS(g, u, distance_back, NULL, BFS_order, NULL, seen)
+    D = distance_back[u]
+
+    for v in range(n):
+        aux_max = 0
+        aux_min = UINT32_MAX
+        for comp in range(nscc):
+            for idx in range(g.neighbors[v]-g.neighbors[v+1]):
+                u = g.neighbors[v][idx]
+                if scc[v] == scc[u]:
+                    aux_min = min(ecc[u] + 1, aux_min)
+            aux_max = max(aux_max, aux_min)
+
+        ecc[v] = min(ecc[v], aux_max)
+
+        if ecc[v] <= D:
+            continue
+
+        # SearchAndBound(G, v)
+        simple_BFS(g, v, distance_from, NULL, BFS_order, NULL, seen)
+
+        aux_ecc = 0
+        u = bitset_first(seen)
+        while u != UINT32_MAX:
+            aux_ecc = max(distance_from[u], aux_ecc)
+            u = bitset_next(seen, u+1)
+
+        ecc[v] = aux_ecc
+        D = max(D, aux_ecc)
+
+        simple_BFS(rev_g, v, distance_back, NULL, BFS_order, scc, seen) # Different graph
+        u = bitset_first(seen)
+        while u != UINT32_MAX:
+            ecc[u] = min(ecc[u], distance_back[u] + aux_ecc)
+            u = bitset_next(seen, u+1)
+
+    return D
+
+def TYY_diameter(G):
+    """
+    TESTS:
+        sage: from sage.graphs.base.static_sparse_graph import TYY_diameter
+        sage: G = Graph([[1,2], [2,3], [3,4], [4,2]])
+        sage: assert(TYY_diameter(G) == 2)
+        sage: G = Graph([[1,2], [2,3], [3,4], [4,1]])
+        sage: assert(TYY_diameter(G) == 2)
+        sage: G = Graph([[1,2], [2,3], [3,4], [4,5], [5,1]])
+        sage: assert(TYY_diameter(G) == 2)
+        sage: G = DiGraph([[1,1], [2,2], [3,3], [4,4]], loops=True)
+        sage: assert(TYY_diameter(G) == 0)
+        sage: G = DiGraph([[1,2], [2,3], [3,1], [1,3], [2,4]])
+        sage: assert(TYY_diameter(G) == 3)
+    """
+    from sage.graphs.graph import Graph 
+    from sage.graphs.digraph import DiGraph 
+
+    if not (isinstance(G, Graph) or isinstance(G, DiGraph)):
+        raise ValueError("G must be either a graph or a digraph")
+
+    cdef MemoryAllocator mem = MemoryAllocator()
+    cdef list int_to_vertex = list(G)
+    cdef short_digraph g
+    init_short_digraph(g, G, edge_labelled=False, vertex_list=int_to_vertex)
+    cdef int *ecc = <int*> mem.malloc(g.n * sizeof(int))
+    sig_on()
+    cdef int D = TYY_diameter_C(g)
+    sig_off()
+    free_short_digraph(g)
+
+    return D
+
 cdef void strongly_connected_components_digraph_C(short_digraph g, int nscc, int *scc, short_digraph output):
     r"""
     Compute the strongly connected components (SCCs) digraph of `g`.

src/sage/graphs/centrality.pyx

@@ -943,7 +943,7 @@ def centrality_closeness_random_k(G, int k=1):
         # Run BFS for random k vertices
         for i in range(k):
             farness = 0
-            simple_BFS(sd, l[i], distance, NULL, waiting_list, seen)
+            simple_BFS(sd, l[i], distance, NULL, waiting_list, NULL, seen)
             for j in range(n):
                 farness += distance[j]
                 partial_farness[j] += distance[j]

src/sage/graphs/digraph.py

@@ -3674,6 +3674,7 @@ def _singleton_in_branching():
     # Aliases to functions defined in other modules
     from sage.graphs.comparability import is_transitive
     from sage.graphs.base.static_sparse_graph import tarjan_strongly_connected_components as strongly_connected_components
+    from sage.graphs.base.static_sparse_graph import TYY_diameter as diameter 
     from sage.graphs.connectivity import is_strongly_connected
     from sage.graphs.connectivity import strongly_connected_components_digraph
     from sage.graphs.connectivity import strongly_connected_components_subgraphs

src/sage/graphs/distances_all_pairs.pyx

@@ -772,7 +772,7 @@ cdef uint32_t * c_eccentricity_bounding(G, vertex_list=None) except NULL:
         cpt += 1
 
         # Compute the exact eccentricity of v
-        LB[v] = simple_BFS(sd, v, distances, NULL, waiting_list, seen)
+        LB[v] = simple_BFS(sd, v, distances, NULL, waiting_list, NULL, seen)
 
         if LB[v] == UINT32_MAX:
             # The graph is not connected. We set maximum value and exit.
@@ -968,7 +968,7 @@ cdef uint32_t diameter_lower_bound_2sweep(short_digraph g,
     cdef uint32_t LB, i, k, tmp
 
     # We do a first BFS from source and get the eccentricity of source
-    LB = simple_BFS(g, source, distances, NULL, waiting_list, seen)
+    LB = simple_BFS(g, source, distances, NULL, waiting_list, NULL, seen)
 
     # If the eccentricity of the source is infinite (very large number), the
     # graph is not connected and so its diameter is infinite.
@@ -977,7 +977,7 @@ cdef uint32_t diameter_lower_bound_2sweep(short_digraph g,
 
     # Then we perform a second BFS from the last visited vertex
     source = waiting_list[g.n - 1]
-    LB = simple_BFS(g, source, distances, predecessors, waiting_list, seen)
+    LB = simple_BFS(g, source, distances, predecessors, waiting_list, NULL, seen)
 
     # We return the computed lower bound
     return LB
@@ -1215,7 +1215,7 @@ cdef uint32_t diameter_iFUB(short_digraph g,
     # We order the vertices by decreasing layers. This is the inverse order of a
     # BFS from m, and so the inverse order of array waiting_list. Distances are
     # stored in array layer.
-    LB = simple_BFS(g, m, layer, NULL, waiting_list, seen)
+    LB = simple_BFS(g, m, layer, NULL, waiting_list, NULL, seen)
     for i in range(n):
         order[i] = waiting_list[n - i - 1]
 
@@ -1241,7 +1241,7 @@ cdef uint32_t diameter_iFUB(short_digraph g,
     # eccentricity already found.
     i = 0
     while 2 * layer[order[i]] > LB and i < n:
-        tmp = simple_BFS(g, order[i], distances, NULL, waiting_list, seen)
+        tmp = simple_BFS(g, order[i], distances, NULL, waiting_list, NULL, seen)
         i += 1
 
         # We update the lower bound

src/sage/graphs/graph.py

@@ -8670,6 +8670,7 @@ def arboricity(self, certificate=False):
     from sage.graphs.graph_decompositions.graph_products import is_cartesian_product
     from sage.graphs.distances_all_pairs import is_distance_regular
     from sage.graphs.base.static_dense_graph import is_strongly_regular
+    from sage.graphs.base.static_sparse_graph import TYY_diameter as diameter
     from sage.graphs.line_graph import is_line_graph
     from sage.graphs.tutte_polynomial import tutte_polynomial
     from sage.graphs.lovasz_theta import lovasz_theta