Ticket #11754: trac_11754.patch

doc/en/reference/graphs.rst

@@ -51 +51 @@
    sage/graphs/pq_trees
    sage/graphs/matchpoly
    sage/graphs/graph_decompositions/vertex_separation
+   sage/graphs/graph_decompositions/rankwidth
+   sage/graphs/modular_decomposition/modular_decomposition
    sage/graphs/convexity_properties
    sage/graphs/distances_all_pairs
    sage/graphs/graph_latex

module_list.py

@@ -405 +405 @@
                          'sage/graphs/planarity/platformTime.h',
                          'sage/graphs/planarity/stack.h']),
 
+    Extension('sage.graphs.graph_decompositions.rankwidth',
+              sources = ['sage/graphs/graph_decompositions/rankwidth.pyx']),
+
     Extension('sage.graphs.spanning_tree',
               sources = ['sage/graphs/spanning_tree.pyx']),
 

sage/graphs/all.py

@@ -12 +12 @@
 import sage.graphs.graph
 import sage.graphs.digraph
 import graph_coloring
+import sage.graphs.graph_decompositions
+import sage.graphs.modular_decomposition.modular_decomposition
 from sage.graphs.cliquer import *
 from graph_database import graph_db_info
 from graph_editor import graph_editor

sage/graphs/graph.py

@@ -2917 +2917 @@
 
         return rs_dict
 
+
+    def rank_decomposition(self, verbose = False):
+        """
+        Returns an rank-decomposition of ``self`` achieving optiml rank-width.
+
+        See the documentation of the ``rankwidth`` module
+        :class:`<rankwidth sage.graphs.graph_decompositions.rankwidth>`.
+
+        INPUT:
+
+        - ``verbose`` (boolean) -- whether to display progress information while
+        computing the decomposition.
+
+        OUTPUT:
+
+        A pair ``(rankwidth, decomposition_tree)``, where ``rankwidth`` is a
+        numerical value and ``decomposition_tree`` is a ternary tree describing
+        the decomposition (cf. the module's documentation).
+
+        See the documentation of the ``rankwidth`` module for more information
+        on the tree
+        :class:`rankwidth <sage.graphs.graph_decompositions.rankwidth.`.
+
+        .. WARNING::
+
+            The current implementation cannot handle graphs of more than 32 vertices.
+
+        EXAMPLE::
+
+            sage: g = graphs.PetersenGraph()
+            sage: rw, tree = g.rank_decomposition()
+            sage: rw
+            3
+            sage: tree
+            Graph on 19 vertices
+            sage: tree.is_tree()
+            True
+        """
+
+        from sage.graphs.graph_decompositions.rankwidth import rank_decomposition
+        return rank_decomposition(self, verbose = verbose)
+
     ### Matching
 
     def matching_polynomial(self, complement=True, name=None):

sage/graphs/graph_decompositions/__init__.py

@@ -1 +1 @@
 # This file is not empty !
 import sage.graphs.graph_decompositions.vertex_separation
+import sage.graphs.graph_decompositions.rankwidth

new file sage/graphs/graph_decompositions/rankwidth.pxd

@@ +1 @@
+cdef extern from *:
+    ctypedef int uint_fast8_t "uint_fast8_t"
+    ctypedef int uint_fast8 "uint_fast8"
+    ctypedef int uint_fast32_t "uint_fast32_t"
+    ctypedef int subset_t "uint_least32_t"
+
+cdef extern from "rankwidth/rw.h":
+    int init_rw_dec(uint_fast8_t n)
+    void destroy_rw()
+    void calculate_level(uint_fast8_t subset_size)
+    uint_fast8_t get_rw()
+    subset_t *cslots
+    subset_t *slots
+    subset_t *adjacency_matrix
+
+cdef extern from "rankwidth/rw.c":
+    uint_fast8_t subset_size
+    uint_fast8_t num_vertices
+
+
+cdef void set_am(int i, int j, int val)
+cdef void print_rank_dec(subset_t s, int l)
+
+
+

new file sage/graphs/graph_decompositions/rankwidth.pyx

@@ +1 @@
+r"""
+Rank Decompositions of graphs.
+
+This modules wraps a C code from Philipp Klaus Krause computing a an optimal
+rank-decomposition [RWKlause]_.
+
+**Definitions :**
+
+Given a graph `G` and a subset `S\subseteq V(G)` of vertices, the *rank-width*
+of `S` in `G`, denoted `rw_G(S)`, is equal to the rank in `GF(2)` of the `|S|
+\times (|V|-|S|)` matrix of the adjacencies between the vertices of `S` and
+`V\backslash S`. By definition, `rw_G(S)` is qual to `rw_G(\overline S)` where
+`\overline S` is the complement of `S` in `V(G)`.
+
+A *rank-decomposition* of `G` is a tree whose `n` leaves are the elements of
+`V(G)`, and whose internal noes have degree 3. In a tree, ay edge naturally
+corresponds to a bipartition of the vertex set : indeed, the reoal of any edge
+splits the tree into two connected components, thus splitting the set of leaves
+(i.e. vertices of `G`) into two sets. Hence we can define for any edge `e\in
+E(G)` a width equal to the value `rw_G(S)` or `rw_G(\overline S)`, where
+`S,\overline S` is the bipartition obtained from `e`. The *rank-width*
+associated to the whole decomposition is then set to the maximum of the width of
+all the edges it contains.
+
+A *rank-decomposition* is said to be optimal for `G` if it is the decomposition
+achieving the minimal *rank-width*.
+
+**RW -- The original source code :**
+
+RW [RWKlause]_ is a program that calculates rank-width and
+rank-decompositions. It is based on ideas from :
+
+    * "Computing rank-width exactly" by Sang-il Oum [Oum]_
+    * "Sopra una formula numerica" by Ernesto Pascal
+    * "Generation of a Vector from the Lexicographical Index" by B.P. Buckles and M. Lybanon [BL]_ 
+    * "Fast additions on masked integers" by Michael D. Adams and David S. Wise [AW]_
+
+**OUTPUT:**
+
+The rank decomposition is returned as a tree whose vertices are subsets of
+`V(G)`. Its leaves, corresponding to the vertices of `G` are sets of 1 elements,
+i.e. singletons.
+
+::
+
+    sage: g = graphs.PetersenGraph()
+    sage: rw, tree = g.rank_decomposition()
+    sage: all(len(v)==1 for v in tree if tree.degree(v) == 1)
+    True
+
+The internal nodes are sets of the decomposition. This way, it is easy to deduce
+the bipartition associated to an edge from the tree. Indeed, two adjacent
+vertices of the tree are comarable sets : they yield the bipartition obtained
+from the smaller of the two and its complement.
+
+::
+
+    sage: g = graphs.PetersenGraph()
+    sage: rw, tree = g.rank_decomposition()
+    sage: u = Set([8, 9, 3, 7])
+    sage: v = Set([8, 9])
+    sage: tree.has_edge(u,v)
+    True
+    sage: m = min(u,v)
+    sage: bipartition = (m, Set(g.vertices()) - m)
+    sage: bipartition
+    ({8, 9}, {0, 1, 2, 3, 4, 5, 6, 7})
+
+.. WARNING::
+
+    The current implementation cannot handle graphs of `\geq 32` vertices.
+
+EXAMPLE::
+
+        sage: g = graphs.PetersenGraph()
+        sage: g.rank_decomposition()
+        (3, Graph on 19 vertices)
+
+AUTHORS:
+
+- Philipp Klaus Krause : Implementation of the C algorithm [RWKlause]_.
+- Nathann Cohen : Interface with Sage and documentation.
+
+REFERENCES:
+
+  .. [RWKlause] Philipp Klaus Krause -- rw v0.2
+    http://pholia.tdi.informatik.uni-frankfurt.de/~philipp/software/rw.shtml
+
+  .. [Oum] Sang-il Oum
+    Computing rank-width exactly
+    Information Processing Letters, 2008
+    vol. 109, n. 13, p. 745--748
+    Elsevier
+    http://mathsci.kaist.ac.kr/~sangil/pdf/2008exp.pdf
+
+  .. [BL] Buckles, B.P. and Lybanon, M.
+    Algorithm 515: generation of a vector from the lexicographical index
+    ACM Transactions on Mathematical Software (TOMS), 1977
+    vol. 3, n. 2, pages 180--182
+    ACM
+
+  .. [AW] Adams, M.D. and Wise, D.S.
+    Fast additions on masked integers
+    ACM SIGPLAN Notices, 2006
+    vol. 41, n.5, pages 39--45
+    ACM
+    http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.86.1801&rep=rep1&type=pdf
+"""
+
+#*****************************************************************************
+#      Copyright (C) 2011 Nathann Cohen <nathann.cohen@gail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+
+include '../../ext/stdsage.pxi'
+include '../../misc/bitset_pxd.pxi'
+include "../../ext/interrupt.pxi"
+
+cdef list id_to_vertices
+cdef dict vertices_to_id
+
+def rank_decomposition(G, verbose = False):
+    r"""
+    Computes an optiml rank-decomposition of the given graph.
+
+    This function is available as a method of the 
+    :class:`Graph <sage.graphs.graph>` class. See
+    :meth:`rank_decomposition <sage.graphs.graph.rank_decomposition>`.
+
+    INPUT:
+
+    - ``verbose`` (boolean) -- whether to display progress information while
+      computing the decomposition.
+
+    OUTPUT:
+
+    A pair ``(rankwidth, decomposition_tree)``, where ``rankwidth`` is a
+    numerical value and ``decomposition_tree`` is a ternary tree describing the
+    decomposition (cf. the module's documentation).
+
+    EXAMPLE::
+
+        sage: from sage.graphs.graph_decompositions.rankwidth import rank_decomposition
+        sage: g = graphs.PetersenGraph()
+        sage: rank_decomposition(g)
+        (3, Graph on 19 vertices)
+
+    On more than 32 vertices::
+
+        sage: g = graphs.RandomGNP(40, .5)
+        sage: rank_decomposition(g)
+        Traceback (most recent call last):
+        ...
+        Exception: The rank decomposition cannot be computed on graphs of >= 32 vertices.
+
+    The empty graph::
+    
+        sage: g = Graph()
+        sage: rank_decomposition(g)
+        (0, Graph on 0 vertices)
+    """
+    cdef int n = G.order()
+
+    if n >= 32:
+        raise Exception("The rank decomposition cannot be computed "+
+                        "on graphs of >= 32 vertices.")
+
+    elif n == 0:
+        from sage.graphs.graph import Graph
+        return (0, Graph())
+
+    global num_vertices
+
+    cdef int i
+
+    if sage_graph_to_matrix(G):
+        raise Exception("There has been a mistake while converting the Sage "+
+                        "graph to a C structure. The memory is probably "+
+                        "insufficient (2^(n+1) is a *LOT*).")
+
+    for 0 <= i < n+1:
+
+        if verbose:
+            print "Calculating for subsets of size ", i, "/",(n+1)
+
+        # We want to properly deal with exceptions, in particular
+        # KeyboardInterrupt. Whatever happens, when this code fails the memory
+        # *HAS* to be freed, as it is particularly greedy (a graph of 29
+        # vertices already requires more than 1GB of memory).
+
+        if not sig_on_no_except():
+            destroy_rw()
+        cython_check_exception()
+
+        # Actual computation
+        calculate_level(i)
+        _sig_off
+
+    cdef int rank_width = <int> get_rw()
+
+    #Original way of displaying the decomposition
+    #print_rank_dec(0x7ffffffful >> (31 - num_vertices), 0)
+    g = mkgraph()
+
+    # Free the memory
+    destroy_rw()
+
+    return (rank_width, g)
+
+cdef int sage_graph_to_matrix(G):
+    r"""
+    Converts the given Sage graph as an adjacency matrix.
+    """
+    global id_to_vertices
+    global vertices_to_id 
+    global adjacency_matrix
+    global slots
+    global cslots
+
+    id_to_vertices = []
+    vertices_to_id = {}
+    
+    global num_vertices
+    num_vertices = G.order()
+
+    cdef int i,j
+
+    # Prepares the C structure for the computation
+    if init_rw_dec(num_vertices):
+        # The original code does not completely frees the memory in case of
+        # error
+        if cslots != NULL:
+            sage_free(cslots)
+        if slots != NULL:
+            sage_free(slots)
+        if adjacency_matrix != NULL:
+            sage_free(adjacency_matrix)
+        return 1
+
+    memset(adjacency_matrix, 0, sizeof(subset_t) * num_vertices)
+
+    # Initializing the lists of vertices
+    for i,v in enumerate(G.vertices()):
+        id_to_vertices.append(v)
+        vertices_to_id[v] = i
+
+    # Filling the matrix
+    for u,v in G.edges(labels = False):
+        if u==v:
+            continue
+        set_am(vertices_to_id[u], vertices_to_id[v], 1)
+
+    # All is fine.
+    return 0
+
+cdef uint_fast32_t bitmask(int i):
+    return(1ul << i)
+
+cdef void set_am(int i, int j, int val):
+    r"""
+    Set/Unset an arc between vertices i and j
+
+    (this function is a copy of what can be found in rankwidth/rw.c)
+    """
+    global adjacency_matrix
+    
+    adjacency_matrix[i] &= ~bitmask(j)
+    adjacency_matrix[j] &= ~bitmask(i)
+
+    if val:
+        adjacency_matrix[i] |= bitmask(j)
+        adjacency_matrix[j] |= bitmask(i)
+
+cdef void print_rank_dec(subset_t s, int l):
+    r"""
+    Prints the current rank decomposition as a text
+
+    This function is a copy of the C routine printing the rank-decomposition is
+    the original source code. It s not used at the moment, but can still prove
+    useful when debugging the code.
+    """
+    global cslots
+
+    print ('\t'*l),
+
+    print "cslot: ", <unsigned int> s
+    if cslots[s] == 0:
+        return
+    print_rank_dec(cslots[s], l + 1)
+    print_rank_dec(s & ~cslots[s], l + 1)
+
+def mkgraph():
+    r"""
+    Returns the graph corresponding the the current rank-decomposition.
+
+    (This function is for interna use)
+
+    EXAMPLE::
+
+        sage: from sage.graphs.graph_decompositions.rankwidth import rank_decomposition
+        sage: g = graphs.PetersenGraph()
+        sage: rank_decomposition(g)
+        (3, Graph on 19 vertices)
+    """
+    global cslots
+    global num_vertices
+    global id_to_vertices
+
+    cdef subset_t s
+    
+    from sage.graphs.graph import Graph
+    g = Graph()
+
+    cdef subset_t * tab = <subset_t *> sage_malloc(sizeof(subset_t) * (2*num_vertices -1))
+    tab[0] = 0x7ffffffful >> (31 - num_vertices)
+
+    cdef int beg = 0
+    cdef int end = 1
+
+    while beg != end:
+
+        s = tab[beg]
+        beg += 1
+
+        if cslots[s] == 0:
+            continue
+
+        g.add_edge(bitset_to_vertex_set(s), bitset_to_vertex_set(s&~cslots[s]))
+        g.add_edge(bitset_to_vertex_set(s), bitset_to_vertex_set(cslots[s]))
+
+        tab[end] = s&~cslots[s]
+        end += 1
+        tab[end] = cslots[s]
+        end += 1
+
+    sage_free(tab)
+    return g
+
+cdef bitset_to_vertex_set(subset_t s):
+    """
+    Returns as a Set object the set corresponding to the given subset_t
+    variable.
+    """
+    from sage.rings.integer import Integer
+    from sage.sets.set import Set
+    from sage.misc.bitset import FrozenBitset
+    return Set(list(FrozenBitset((Integer(<unsigned int> s).binary())[::-1])))

new file sage/graphs/graph_decompositions/rankwidth/README

@@ +1 @@
+This program calculates rank-width and rank-decompositions. It is based on ideas from "Computing rank-width exactly" by Sang-il Oum, "Sopra una formula numerica" by Ernesto Pascal, "Generation of a Vector from the Lexicographical Index" by B.P. Buckles and M. Lybanon and "Fast additions on masked integers" by Michael D. Adams and David S. Wise.
+On today's computers it works quite well up to graph sizes of about 28 nodes.
+
+Assuming that libigraph is installed it can be compiled using e.g. "gcc -O2 --std=c99 -pedantic rw.c simplerw.c -ligraph".
+
+It supports a number of input graph formats; there is an example graph included, for which rank-width and rank-decomposition can be calculated using "./a.out --edgelist igrid".
+
+rw.h / rw.c provide a general mechanism for calculating rank-width and rank-decompositions. simplerw.c provides a simple command-line interface.
+

new file sage/graphs/graph_decompositions/rankwidth/__init__.py

@@ +1 @@
+# This file is not empty !

new file sage/graphs/graph_decompositions/rankwidth/rw.c

@@ +1 @@
+// Calculate rank-width and rank-decompositions of graphs.
+
+// Philipp Klaus Krause, philipp@informatik.uni-frankfurt.de, pkk@spth.de, 2009 - 2012
+// Copyright (c) 2009-2012 Philipp Klaus Krause
+
+// 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 2 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, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+    
+#include <stdlib.h>
+
+#include "rw.h"
+
+static uint_fast8_t subset_size;
+
+static uint_fast8_t num_vertices;
+
+subset_t *adjacency_matrix;
+
+// This array is meant to contain max{f(Y), w(f, Y)} at the position corresponding to Y.
+static uint_fast8_t *slots;
+subset_t *cslots;
+
+uint_fast8_t cut_rank(const subset_t s)
+{
+        subset_t am[num_vertices];
+        subset_t x, y;
+        uint_fast8_t i, j;
+        uint_fast8_t rank = 0;
+
+        for(i = 0; i < num_vertices; i++)
+                am[i] =  (s & (1ul << i)) ? 0 : (adjacency_matrix[i] & s);
+
+        for(i = 0; i < num_vertices; i++)
+        {
+                y = 0;
+                for(j = rank; j < num_vertices; j++)
+                {
+                        x = am[j];
+                        if(x & 1)
+                        {
+                                if(!y)
+                                {
+                                        y = x;
+                                        am[j] = am[rank++];
+                                        continue;
+                                }
+                                x ^= y;
+                        }
+                        am[j] = x >> 1;
+                }
+        }
+
+        return(rank);
+}
+
+// Compute the binomial coefficient C(n, k).
+// Time complexity: O(n) * complexity of integer division.
+// Could be replaced by a lookup table of size O(n^2), but since this is not the bopttleneck, it doesn't matter.
+subset_t binomial_coefficient(uint_fast8_t n, uint_fast8_t k)
+{
+        uint_fast8_t i, delta, max;
+        subset_t c;
+
+        if(n < k)
+                return(0);
+
+        if(n == k)
+                return(1);
+
+        if(k < n - k)
+        {
+                delta = n - k;
+                max = k;
+        }
+        else
+        {
+                delta = k;
+                max = n - k;
+        }
+
+        c = delta + 1;
+
+        for(i = 2; i <= max; i++)
+                c = (c * (delta + i)) / i;
+
+        return(c);
+}
+
+// Convert unsigned integer from combinadic to machine representation.
+subset_t comb_to_int(subset_t c)
+{
+        uint_fast8_t k, n;
+        subset_t i;
+
+        i = 0;
+        for(k = 0, n = 1; k < num_vertices; k++, c >>= 1)
+                if(c & 1)
+                        i += binomial_coefficient(k, n++);
+        return(i);
+}
+
+// Return largest value v where v < a and  Choose(v,b) <= x.
+static uint_fast8_t largest_v(uint_fast8_t a, uint_fast8_t b, subset_t x)
+{
+        uint_fast8_t v = a - 1;
+           
+        while(binomial_coefficient(v,b) > x)
+                v--;
+
+        return(v);
+}
+
+// Convert unsigned integer from machine representation to combinadic.
+static subset_t int_to_comb(subset_t i)
+{
+        uint_fast8_t j, v;
+        uint_fast8_t a = num_vertices;
+        uint_fast8_t b = subset_size;
+        subset_t c = 0;
+
+        for(j = 0; j < subset_size; j++, b--)
+        {
+                v = largest_v(a, b, i);
+                i = i - binomial_coefficient(v, b);
+                a = v;
+                c |= (1ul << v);
+        }
+
+        return(c);
+}
+
+// Masked increment.
+static subset_t subset_inc(subset_t v, subset_t mask)
+{
+        return((v - mask) & mask);
+}
+
+// Returns rank-width for a subset of size at least 2 given that slots already contains correct values for all nonempty subsets of sizes less than the size of s.
+// This is where most of the time is spent.
+uint_fast8_t width(subset_t s)
+{
+        uint_fast8_t w = UINT_FAST8_MAX, v, v1, v2;
+        subset_t ss;
+        subset_t cs;
+        for(ss = subset_inc(0, s); ss != s; ss = subset_inc(ss, s))
+        {
+                v1 = slots[ss];
+                v2 = slots[s & ~ss];
+                v = v1 > v2 ? v1 : v2;
+                if(v < w)
+                {
+                        w = v;
+                        cs = ss;
+                }
+        }
+        cslots[s] = cs;
+        return(w);
+}
+
+void fill_slot(subset_t i)
+{
+        uint_fast8_t v, w;
+        subset_t s = int_to_comb(i);
+        v = cut_rank(s);
+        w = width(s);
+        slots[s] = v > w ? v : w;
+}
+
+void calculate_level(uint_fast8_t ss)
+{
+        uint_fast8_t i;
+
+        subset_size = ss;
+
+        if(subset_size == 0)
+                slots[0] = 0;
+        else if(subset_size == 1)
+                for(i = 0; i < num_vertices; i++)
+                {
+                        slots[1ul << i] = cut_rank(1ul << i);
+                        cslots[1ul << i] = 0x0;
+                }
+        else
+        {
+                subset_t i;
+                const subset_t end = binomial_coefficient(num_vertices, subset_size);
+                for(i = 0; i < end; i++)
+                        fill_slot(i);
+        }
+}
+
+void calculate_all(void)
+{
+        uint_fast8_t i;
+
+        for(i = 0; i < num_vertices; i++)
+        {
+                slots[1ul << i] = cut_rank(1ul << i);
+                cslots[1ul << i] = 0x0;
+        }
+
+        for(subset_size = 2; subset_size <= num_vertices; subset_size++)
+        {
+                subset_t i;
+                const subset_t end = binomial_coefficient(num_vertices, subset_size);
+                for(i = 0; i < end; i++)
+                        fill_slot(i);
+        }
+}
+
+int init_rw(uint_fast8_t n)
+{
+        // If sizeof(uint_fast8_t) * (1ul << n) overflows, it wraps around to 0, since size_t and unsigned long are unsigned integer types.
+        if(n > MAX_VERTICES || n && !(sizeof(uint_fast8_t) * (1ul << n)))
+                return(-1);
+
+        num_vertices = n;
+        adjacency_matrix = malloc(sizeof(subset_t) * n);
+        slots = malloc(sizeof(uint_fast8_t) * (1ul << n));
+        cslots = 0;
+        return((adjacency_matrix && slots) ? 0 : -1);
+}
+
+int init_rw_dec(uint_fast8_t n)
+{
+        // If sizeof(subset_t) * (1ul << n) overflows, it wraps around to 0, since size_t and unsigned long are unsigned integer types.
+        if(n && !(sizeof(subset_t) * (1ul << n)))
+                return(-1);
+
+        if(init_rw(n))
+                return(-1);
+        cslots = malloc(sizeof(subset_t) * (1ul << n));
+        return(cslots ? 0 : -1);
+}
+
+void destroy_rw(void)
+{
+        free(cslots);
+        free(slots);
+        free(adjacency_matrix);
+}
+
+uint_fast8_t get_rw(void)
+{
+        return(slots[0xfffffffful >> (32 - num_vertices)]);
+}
+

new file sage/graphs/graph_decompositions/rankwidth/rw.h

@@ +1 @@
+// Calculate rank-width and rank-decompositions of graphs.
+
+// Philipp Klaus Krause, philipp@informatik.uni-frankfurt.de, pkk@spth.de, 2009 - 2011
+// Copyright (c) 2009-2011 Philipp Klaus Krause
+
+// 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 2 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, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+// This is a program that calculates rank-width and rank-decompositions. It is based on ideas from "Computing rank-width exactly" by Sang-il Oum, "Sopra una formula numerica" by Ernesto Pascal and "Generation of a Vector from the Lexicographical Index" by B.P. Buckles and M. Lybanon.
+// On 2009's computers it works quite well up to graph sizes of about 28 nodes. 
+// It is an implementation of the trivial algorithm from "Computing rank-width exactly". For larger graphs (more than about 40 nodes) the more algorithm based on fast subset convolution from the same paper would probably be faster.
+
+#include <stdint.h>
+
+// Use data type uint_leastN_t. N is an upper limit on the size of the graphs that can be handled. N=32 seems to be a good compromise for now (the code works well with other values of N).
+// uint_leastN_t is faster than uint_fastN_t here, since the bottleneck is cache misses.
+typedef uint_least32_t subset_t;
+#define MAX_VERTICES 32
+
+// Input graph.
+//extern subset_t adjacency_matrix[NUM_VERTICES];
+extern subset_t *adjacency_matrix;
+
+// Output rank-decomposition
+extern subset_t *cslots;
+
+// Initialization (for getting rank-width only). Returns 0 on success.
+//int init_rw(uint_fast8_t n);
+
+// Initialization (for getting both rank-width and rank-decomposition). Returns 0 on success.
+int init_rw_dec(uint_fast8_t n);
+
+// Free memory allocated during initialization.
+void destroy_rw(void);
+
+// Calculate everything. May take some time.
+void calculate_all(void);
+
+// Calculate a single level only
+void calculate_level(uint_fast8_t subset_size);
+
+// Get the rank-width.
+uint_fast8_t get_rw(void);
+

new file sage/graphs/graph_decompositions/rankwidth/simplerw.c

@@ +1 @@
+// A simple program for calculating rank-width and rank-decompositions.
+// Compile using gcc -O2 --std=c99 -pedantic rw.c simplerw.c -ligraph
+
+// Philipp Klaus Krause, philipp@informatik.uni-frankfurt.de, pkk@spth.de, 2009 - 2011
+// Copyright (c) 2009-2011 Philipp Klaus Krause
+
+// 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 2 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, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <igraph/igraph.h>
+
+#include "rw.h"
+
+static int num_vertices = 18;
+
+static uint_fast32_t bitmask(int i)
+{
+        return(1ul << i);
+}
+
+static void set_am(int i, int j, int val)
+{
+        adjacency_matrix[i] &= ~bitmask(j);
+        adjacency_matrix[j] &= ~bitmask(i);
+        if(val)
+        {
+                adjacency_matrix[i] |= bitmask(j);
+                adjacency_matrix[j] |= bitmask(i);
+        }
+}
+
+static void tab(int l)
+{
+        while(l--)
+                printf("\t");
+}
+
+void print_rank_dec(subset_t s, int l)
+{
+        tab(l);
+        printf("cslot: %x\n", (unsigned)s);
+        if(cslots[s] == 0)
+                return;
+        print_rank_dec(cslots[s], l + 1);
+        print_rank_dec(s & ~cslots[s], l + 1);
+}
+
+/*int random_graph(void)
+{
+        int i, j;
+
+        if(init_rw_dec(num_vertices))
+                return(-1);
+
+        srand(23);
+        for(i = 0; i < num_vertices; i++)
+                for(j = 0; j < num_vertices; j++)
+                        set_am(i, j, rand() % 2);
+
+        for(i = 0; i < num_vertices; i++)
+                printf("Adj. mat.: %x\n", (unsigned)(adjacency_matrix[i]));
+
+        return(0);
+}*/
+
+int read_graph(const char *format, const char * filename)
+{
+        int ret, i, j;
+        FILE *file;
+        igraph_t igraph;
+        igraph_matrix_t imatrix;
+        igraph_bool_t isimple;
+
+        // Open file
+        if(!(file = fopen(filename, "r")))
+        {
+                printf("Failed to open file %s.\n", filename);
+                return(-1);
+        }
+        
+        // Read graph from file
+        if(!strcmp("--edgelist", format))
+                ret = igraph_read_graph_edgelist(&igraph, file, 0, 0);
+        else if(!strcmp("--graphdb", format))
+                ret = igraph_read_graph_graphdb(&igraph, file, 0);
+        else if(!strcmp("--graphml", format))
+                ret = igraph_read_graph_graphml(&igraph, file, 0);
+        else if(!strcmp("--gml", format))
+                ret = igraph_read_graph_gml(&igraph, file);
+        else if(!strcmp("--pajek", format))
+                ret = igraph_read_graph_pajek(&igraph, file);
+        else
+        {
+                printf("Unknown file format %s.\n", format);
+                fclose(file);
+                return(-1);
+        }
+
+        fclose(file);
+
+        if(ret)
+        {
+                printf("Failed to read a graph from file %s.\n", filename);
+                return(-1);
+        }
+
+        // Check for undirectedness
+        if(igraph_is_directed(&igraph) || igraph_is_simple(&igraph, &isimple) || !isimple)
+        {
+                printf("Input is not a simple, undirected graph from file %s.\n", filename);
+                igraph_destroy(&igraph);
+                return(-1);
+        }
+
+        // Get adjacency matrix
+        if(igraph_matrix_init(&imatrix, 1, 1))
+        {
+                igraph_destroy(&igraph);
+                return(-1);
+        }
+        igraph_get_adjacency(&igraph, &imatrix, IGRAPH_GET_ADJACENCY_BOTH);
+        igraph_destroy(&igraph);
+        if(igraph_matrix_nrow(&imatrix) > MAX_VERTICES)
+        {
+                igraph_matrix_destroy(&imatrix);
+                printf("Graph too large.\n");
+                return(-1);
+        }
+
+        num_vertices = igraph_matrix_nrow(&imatrix);
+        
+        if(init_rw_dec(num_vertices))
+        {
+                destroy_rw();
+                printf("Not enough memory.\n");
+                return(-1);
+        }
+
+        for(i = 0; i < num_vertices; i++)
+                for(j = 0; j < num_vertices; j++)
+                        set_am(i, j, MATRIX(imatrix, i, j));
+
+        igraph_matrix_destroy(&imatrix);
+        
+        return(ret ? -1 : 0);
+}
+
+int main(int argc, char *argv[])
+{
+        int i;
+
+        if(argc /*==*/ <= 2 || argc > 4)
+        {
+                printf("Usage: rw [--format filename]\n");
+                printf("Supported formats: edgelist, graphdb, graphml, gml, pajek.\n");
+                return(-1);
+        }
+
+        if(/*argc <= 1 ? random_graph() :*/ read_graph(argv[1], argv[2]))
+        {
+                printf("Failed to create input graph.\n");
+                return(-1);
+        }
+
+        printf("%d vertices.\n", (int)num_vertices);
+
+        for(i = 0; i <= num_vertices; i++)
+        {
+                printf("Calculating for subsets of size %d.\n", i);
+                calculate_level(i);
+        }
+
+        printf("rank-width: %d\n", (int)get_rw());
+
+        print_rank_dec(0x7ffffffful >> (31 - num_vertices), 0);
+
+        destroy_rw();
+
+        return(0);
+}
+

sage/graphs/modular_decomposition/__init__.py

@@ -1 +1 @@
 # This file is not empty !
+
+import sage.graphs.modular_decomposition.modular_decomposition
+

sage/graphs/modular_decomposition/modular_decomposition.pyx

@@ +1 @@
+r"""
+Modular decomposition
+"""
+
 #####################################################
 # The following code is mainly a Cythonized
 # copy of code found in src/random.c and src/dm.c
@@ -26 +30 @@
 
 cpdef modular_decomposition(g):
     r"""
-    Returns a modular decomposition of the given graph
+    Returns a modular decomposition of the given graph.
 
     INPUT:
 
@@ -42 +46 @@
             * ``"Prime"``
             * ``"Serie"``
 
-        * The list of submodules (as list of pairs ``(type, list)``, recursively...)
-          or the vertex's name if the module is a singleton.          
+        * The list of submodules (as list of pairs ``(type, list)``,
+          recursively...)  or the vertex's name if the module is a
+          singleton.
 
     .. NOTE::
 
-    As this fuction could be used by efficient C routines, the vertices returned are
-    not labels but identifiants from ``[0, ..., g.order()-1]``
-
-
+        As this fuction could be used by efficient C routines, the
+        vertices returned are not labels but identifiants from ``[0,
+        ..., g.order()-1]``
 
     ALGORITHM:
 
     This function uses a C implementation of a 2-step algorithm 
-    implemented by Fabien de Montgolfier [FMDec]_ :
+    implemented by Fabien de Montgolfier [FMDecb]_ :
 
-        * Computation of a factorizing permutation [HabibViennot1999]_.
+        * Computation of a factorizing permutation [HabibViennot1999b]_.
 
-        * Computation of the tree itself [CapHabMont02]_.
-
-    REFERENCE:
-
-    .. [FMDec] Fabien de Montgolfier
-      http://www.liafa.jussieu.fr/~fm/algos/index.html
-
-    .. [HabibViennot1999] Michel Habib, Christiphe Paul, Laurent Viennot
-      Partition refinement techniques: An interesting algorithmic tool kit
-      International Journal of Foundations of Computer Science
-      vol. 10 n2 pp.147--170, 1999
-
-    .. [CapHabMont02] C. Capelle, M. Habib et F. de Montgolfier
-      Graph decomposition and Factorising Permutations
-      Discrete Mathematics and Theoretical Computer Sciences, vol 5 no. 1 , 2002.
+        * Computation of the tree itself [CapHabMont02b]_.
 
     EXAMPLES:
 
@@ -96 +86 @@
         sage: g.add_edge(0,6)
         sage: modular_decomposition(g)
         ('Serie', [0, ('Parallel', [5, ('Serie', [1, 4, 3, 2]), 6])])
+
+
+    REFERENCES:
+
+    .. [FMDecb] Fabien de Montgolfier
+      http://www.liafa.jussieu.fr/~fm/algos/index.html
+
+    .. [HabibViennot1999b] Michel Habib, Christiphe Paul, Laurent Viennot
+      Partition refinement techniques: An interesting algorithmic tool kit
+      International Journal of Foundations of Computer Science
+      vol. 10 n2 pp.147--170, 1999
+
+    .. [CapHabMont02b] C. Capelle, M. Habib et F. de Montgolfier
+      Graph decomposition and Factorising Permutations
+      Discrete Mathematics and Theoretical Computer Sciences, vol 5 no. 1 , 2002.
     """
     cdef c_graphe G
     cdef c_adj * a

setup.py

@@ -905 +905 @@
                      'sage.graphs.base',
                      'sage.graphs.modular_decomposition',
                      'sage.graphs.graph_decompositions',
+                     'sage.graphs.graph_decompositions',
                      
                      'sage.groups',
                      'sage.groups.abelian_gps',