Ticket #11754: trac_11754.patch

File trac_11754.patch, 33.7 KB (added by ncohen, 9 years ago)
  • doc/en/reference/graphs.rst

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1314471808 -7200
    # Node ID a494ee165a85263c4db159eea432362cfc9b8f8e
    # Parent  22974c396744c772d86d0496aa98a91d6bc09ccf
    trac 11754 -- Computation of rank-decompositions in Sage
    
    diff --git a/doc/en/reference/graphs.rst b/doc/en/reference/graphs.rst
    a b  
    5151   sage/graphs/pq_trees
    5252   sage/graphs/matchpoly
    5353   sage/graphs/graph_decompositions/vertex_separation
     54   sage/graphs/graph_decompositions/rankwidth
     55   sage/graphs/modular_decomposition/modular_decomposition
    5456   sage/graphs/convexity_properties
    5557   sage/graphs/distances_all_pairs
    5658   sage/graphs/graph_latex
  • module_list.py

    diff --git a/module_list.py b/module_list.py
    a b  
    405405                         'sage/graphs/planarity/platformTime.h',
    406406                         'sage/graphs/planarity/stack.h']),
    407407
     408    Extension('sage.graphs.graph_decompositions.rankwidth',
     409              sources = ['sage/graphs/graph_decompositions/rankwidth.pyx']),
     410
    408411    Extension('sage.graphs.spanning_tree',
    409412              sources = ['sage/graphs/spanning_tree.pyx']),
    410413
  • sage/graphs/all.py

    diff --git a/sage/graphs/all.py b/sage/graphs/all.py
    a b  
    1212import sage.graphs.graph
    1313import sage.graphs.digraph
    1414import graph_coloring
     15import sage.graphs.graph_decompositions
     16import sage.graphs.modular_decomposition.modular_decomposition
    1517from sage.graphs.cliquer import *
    1618from graph_database import graph_db_info
    1719from graph_editor import graph_editor
  • sage/graphs/graph.py

    diff --git a/sage/graphs/graph.py b/sage/graphs/graph.py
    a b  
    29172917
    29182918        return rs_dict
    29192919
     2920
     2921    def rank_decomposition(self, verbose = False):
     2922        """
     2923        Returns an rank-decomposition of ``self`` achieving optiml rank-width.
     2924
     2925        See the documentation of the ``rankwidth`` module
     2926        :class:`<rankwidth sage.graphs.graph_decompositions.rankwidth>`.
     2927
     2928        INPUT:
     2929
     2930        - ``verbose`` (boolean) -- whether to display progress information while
     2931        computing the decomposition.
     2932
     2933        OUTPUT:
     2934
     2935        A pair ``(rankwidth, decomposition_tree)``, where ``rankwidth`` is a
     2936        numerical value and ``decomposition_tree`` is a ternary tree describing
     2937        the decomposition (cf. the module's documentation).
     2938
     2939        See the documentation of the ``rankwidth`` module for more information
     2940        on the tree
     2941        :class:`rankwidth <sage.graphs.graph_decompositions.rankwidth.`.
     2942
     2943        .. WARNING::
     2944
     2945            The current implementation cannot handle graphs of more than 32 vertices.
     2946
     2947        EXAMPLE::
     2948
     2949            sage: g = graphs.PetersenGraph()
     2950            sage: rw, tree = g.rank_decomposition()
     2951            sage: rw
     2952            3
     2953            sage: tree
     2954            Graph on 19 vertices
     2955            sage: tree.is_tree()
     2956            True
     2957        """
     2958
     2959        from sage.graphs.graph_decompositions.rankwidth import rank_decomposition
     2960        return rank_decomposition(self, verbose = verbose)
     2961
    29202962    ### Matching
    29212963
    29222964    def matching_polynomial(self, complement=True, name=None):
  • sage/graphs/graph_decompositions/__init__.py

    diff --git a/sage/graphs/graph_decompositions/__init__.py b/sage/graphs/graph_decompositions/__init__.py
    a b  
    11# This file is not empty !
    22import sage.graphs.graph_decompositions.vertex_separation
     3import sage.graphs.graph_decompositions.rankwidth
  • new file sage/graphs/graph_decompositions/rankwidth.pxd

    diff --git a/sage/graphs/graph_decompositions/rankwidth.pxd b/sage/graphs/graph_decompositions/rankwidth.pxd
    new file mode 100644
    - +  
     1cdef extern from *:
     2    ctypedef int uint_fast8_t "uint_fast8_t"
     3    ctypedef int uint_fast8 "uint_fast8"
     4    ctypedef int uint_fast32_t "uint_fast32_t"
     5    ctypedef int subset_t "uint_least32_t"
     6
     7cdef extern from "rankwidth/rw.h":
     8    int init_rw_dec(uint_fast8_t n)
     9    void destroy_rw()
     10    void calculate_level(uint_fast8_t subset_size)
     11    uint_fast8_t get_rw()
     12    subset_t *cslots
     13    subset_t *slots
     14    subset_t *adjacency_matrix
     15
     16cdef extern from "rankwidth/rw.c":
     17    uint_fast8_t subset_size
     18    uint_fast8_t num_vertices
     19
     20
     21cdef void set_am(int i, int j, int val)
     22cdef void print_rank_dec(subset_t s, int l)
     23
     24
     25
  • new file sage/graphs/graph_decompositions/rankwidth.pyx

    diff --git a/sage/graphs/graph_decompositions/rankwidth.pyx b/sage/graphs/graph_decompositions/rankwidth.pyx
    new file mode 100644
    - +  
     1r"""
     2Rank Decompositions of graphs.
     3
     4This modules wraps a C code from Philipp Klaus Krause computing a an optimal
     5rank-decomposition [RWKlause]_.
     6
     7**Definitions :**
     8
     9Given a graph `G` and a subset `S\subseteq V(G)` of vertices, the *rank-width*
     10of `S` in `G`, denoted `rw_G(S)`, is equal to the rank in `GF(2)` of the `|S|
     11\times (|V|-|S|)` matrix of the adjacencies between the vertices of `S` and
     12`V\backslash S`. By definition, `rw_G(S)` is qual to `rw_G(\overline S)` where
     13`\overline S` is the complement of `S` in `V(G)`.
     14
     15A *rank-decomposition* of `G` is a tree whose `n` leaves are the elements of
     16`V(G)`, and whose internal noes have degree 3. In a tree, ay edge naturally
     17corresponds to a bipartition of the vertex set : indeed, the reoal of any edge
     18splits the tree into two connected components, thus splitting the set of leaves
     19(i.e. vertices of `G`) into two sets. Hence we can define for any edge `e\in
     20E(G)` a width equal to the value `rw_G(S)` or `rw_G(\overline S)`, where
     21`S,\overline S` is the bipartition obtained from `e`. The *rank-width*
     22associated to the whole decomposition is then set to the maximum of the width of
     23all the edges it contains.
     24
     25A *rank-decomposition* is said to be optimal for `G` if it is the decomposition
     26achieving the minimal *rank-width*.
     27
     28**RW -- The original source code :**
     29
     30RW [RWKlause]_ is a program that calculates rank-width and
     31rank-decompositions. It is based on ideas from :
     32
     33    * "Computing rank-width exactly" by Sang-il Oum [Oum]_
     34    * "Sopra una formula numerica" by Ernesto Pascal
     35    * "Generation of a Vector from the Lexicographical Index" by B.P. Buckles and M. Lybanon [BL]_
     36    * "Fast additions on masked integers" by Michael D. Adams and David S. Wise [AW]_
     37
     38**OUTPUT:**
     39
     40The rank decomposition is returned as a tree whose vertices are subsets of
     41`V(G)`. Its leaves, corresponding to the vertices of `G` are sets of 1 elements,
     42i.e. singletons.
     43
     44::
     45
     46    sage: g = graphs.PetersenGraph()
     47    sage: rw, tree = g.rank_decomposition()
     48    sage: all(len(v)==1 for v in tree if tree.degree(v) == 1)
     49    True
     50
     51The internal nodes are sets of the decomposition. This way, it is easy to deduce
     52the bipartition associated to an edge from the tree. Indeed, two adjacent
     53vertices of the tree are comarable sets : they yield the bipartition obtained
     54from the smaller of the two and its complement.
     55
     56::
     57
     58    sage: g = graphs.PetersenGraph()
     59    sage: rw, tree = g.rank_decomposition()
     60    sage: u = Set([8, 9, 3, 7])
     61    sage: v = Set([8, 9])
     62    sage: tree.has_edge(u,v)
     63    True
     64    sage: m = min(u,v)
     65    sage: bipartition = (m, Set(g.vertices()) - m)
     66    sage: bipartition
     67    ({8, 9}, {0, 1, 2, 3, 4, 5, 6, 7})
     68
     69.. WARNING::
     70
     71    The current implementation cannot handle graphs of `\geq 32` vertices.
     72
     73EXAMPLE::
     74
     75        sage: g = graphs.PetersenGraph()
     76        sage: g.rank_decomposition()
     77        (3, Graph on 19 vertices)
     78
     79AUTHORS:
     80
     81- Philipp Klaus Krause : Implementation of the C algorithm [RWKlause]_.
     82- Nathann Cohen : Interface with Sage and documentation.
     83
     84REFERENCES:
     85
     86  .. [RWKlause] Philipp Klaus Krause -- rw v0.2
     87    http://pholia.tdi.informatik.uni-frankfurt.de/~philipp/software/rw.shtml
     88
     89  .. [Oum] Sang-il Oum
     90    Computing rank-width exactly
     91    Information Processing Letters, 2008
     92    vol. 109, n. 13, p. 745--748
     93    Elsevier
     94    http://mathsci.kaist.ac.kr/~sangil/pdf/2008exp.pdf
     95
     96  .. [BL] Buckles, B.P. and Lybanon, M.
     97    Algorithm 515: generation of a vector from the lexicographical index
     98    ACM Transactions on Mathematical Software (TOMS), 1977
     99    vol. 3, n. 2, pages 180--182
     100    ACM
     101
     102  .. [AW] Adams, M.D. and Wise, D.S.
     103    Fast additions on masked integers
     104    ACM SIGPLAN Notices, 2006
     105    vol. 41, n.5, pages 39--45
     106    ACM
     107    http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.86.1801&rep=rep1&type=pdf
     108"""
     109
     110#*****************************************************************************
     111#      Copyright (C) 2011 Nathann Cohen <nathann.cohen@gail.com>
     112#
     113# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     114#                         http://www.gnu.org/licenses/
     115#*****************************************************************************
     116
     117
     118include '../../ext/stdsage.pxi'
     119include '../../misc/bitset_pxd.pxi'
     120include "../../ext/interrupt.pxi"
     121
     122cdef list id_to_vertices
     123cdef dict vertices_to_id
     124
     125def rank_decomposition(G, verbose = False):
     126    r"""
     127    Computes an optiml rank-decomposition of the given graph.
     128
     129    This function is available as a method of the
     130    :class:`Graph <sage.graphs.graph>` class. See
     131    :meth:`rank_decomposition <sage.graphs.graph.rank_decomposition>`.
     132
     133    INPUT:
     134
     135    - ``verbose`` (boolean) -- whether to display progress information while
     136      computing the decomposition.
     137
     138    OUTPUT:
     139
     140    A pair ``(rankwidth, decomposition_tree)``, where ``rankwidth`` is a
     141    numerical value and ``decomposition_tree`` is a ternary tree describing the
     142    decomposition (cf. the module's documentation).
     143
     144    EXAMPLE::
     145
     146        sage: from sage.graphs.graph_decompositions.rankwidth import rank_decomposition
     147        sage: g = graphs.PetersenGraph()
     148        sage: rank_decomposition(g)
     149        (3, Graph on 19 vertices)
     150
     151    On more than 32 vertices::
     152
     153        sage: g = graphs.RandomGNP(40, .5)
     154        sage: rank_decomposition(g)
     155        Traceback (most recent call last):
     156        ...
     157        Exception: The rank decomposition cannot be computed on graphs of >= 32 vertices.
     158
     159    The empty graph::
     160   
     161        sage: g = Graph()
     162        sage: rank_decomposition(g)
     163        (0, Graph on 0 vertices)
     164    """
     165    cdef int n = G.order()
     166
     167    if n >= 32:
     168        raise Exception("The rank decomposition cannot be computed "+
     169                        "on graphs of >= 32 vertices.")
     170
     171    elif n == 0:
     172        from sage.graphs.graph import Graph
     173        return (0, Graph())
     174
     175    global num_vertices
     176
     177    cdef int i
     178
     179    if sage_graph_to_matrix(G):
     180        raise Exception("There has been a mistake while converting the Sage "+
     181                        "graph to a C structure. The memory is probably "+
     182                        "insufficient (2^(n+1) is a *LOT*).")
     183
     184    for 0 <= i < n+1:
     185
     186        if verbose:
     187            print "Calculating for subsets of size ", i, "/",(n+1)
     188
     189        # We want to properly deal with exceptions, in particular
     190        # KeyboardInterrupt. Whatever happens, when this code fails the memory
     191        # *HAS* to be freed, as it is particularly greedy (a graph of 29
     192        # vertices already requires more than 1GB of memory).
     193
     194        if not sig_on_no_except():
     195            destroy_rw()
     196        cython_check_exception()
     197
     198        # Actual computation
     199        calculate_level(i)
     200        _sig_off
     201
     202    cdef int rank_width = <int> get_rw()
     203
     204    #Original way of displaying the decomposition
     205    #print_rank_dec(0x7ffffffful >> (31 - num_vertices), 0)
     206    g = mkgraph()
     207
     208    # Free the memory
     209    destroy_rw()
     210
     211    return (rank_width, g)
     212
     213cdef int sage_graph_to_matrix(G):
     214    r"""
     215    Converts the given Sage graph as an adjacency matrix.
     216    """
     217    global id_to_vertices
     218    global vertices_to_id
     219    global adjacency_matrix
     220    global slots
     221    global cslots
     222
     223    id_to_vertices = []
     224    vertices_to_id = {}
     225   
     226    global num_vertices
     227    num_vertices = G.order()
     228
     229    cdef int i,j
     230
     231    # Prepares the C structure for the computation
     232    if init_rw_dec(num_vertices):
     233        # The original code does not completely frees the memory in case of
     234        # error
     235        if cslots != NULL:
     236            sage_free(cslots)
     237        if slots != NULL:
     238            sage_free(slots)
     239        if adjacency_matrix != NULL:
     240            sage_free(adjacency_matrix)
     241        return 1
     242
     243    memset(adjacency_matrix, 0, sizeof(subset_t) * num_vertices)
     244
     245    # Initializing the lists of vertices
     246    for i,v in enumerate(G.vertices()):
     247        id_to_vertices.append(v)
     248        vertices_to_id[v] = i
     249
     250    # Filling the matrix
     251    for u,v in G.edges(labels = False):
     252        if u==v:
     253            continue
     254        set_am(vertices_to_id[u], vertices_to_id[v], 1)
     255
     256    # All is fine.
     257    return 0
     258
     259cdef uint_fast32_t bitmask(int i):
     260    return(1ul << i)
     261
     262cdef void set_am(int i, int j, int val):
     263    r"""
     264    Set/Unset an arc between vertices i and j
     265
     266    (this function is a copy of what can be found in rankwidth/rw.c)
     267    """
     268    global adjacency_matrix
     269   
     270    adjacency_matrix[i] &= ~bitmask(j)
     271    adjacency_matrix[j] &= ~bitmask(i)
     272
     273    if val:
     274        adjacency_matrix[i] |= bitmask(j)
     275        adjacency_matrix[j] |= bitmask(i)
     276
     277cdef void print_rank_dec(subset_t s, int l):
     278    r"""
     279    Prints the current rank decomposition as a text
     280
     281    This function is a copy of the C routine printing the rank-decomposition is
     282    the original source code. It s not used at the moment, but can still prove
     283    useful when debugging the code.
     284    """
     285    global cslots
     286
     287    print ('\t'*l),
     288
     289    print "cslot: ", <unsigned int> s
     290    if cslots[s] == 0:
     291        return
     292    print_rank_dec(cslots[s], l + 1)
     293    print_rank_dec(s & ~cslots[s], l + 1)
     294
     295def mkgraph():
     296    r"""
     297    Returns the graph corresponding the the current rank-decomposition.
     298
     299    (This function is for interna use)
     300
     301    EXAMPLE::
     302
     303        sage: from sage.graphs.graph_decompositions.rankwidth import rank_decomposition
     304        sage: g = graphs.PetersenGraph()
     305        sage: rank_decomposition(g)
     306        (3, Graph on 19 vertices)
     307    """
     308    global cslots
     309    global num_vertices
     310    global id_to_vertices
     311
     312    cdef subset_t s
     313   
     314    from sage.graphs.graph import Graph
     315    g = Graph()
     316
     317    cdef subset_t * tab = <subset_t *> sage_malloc(sizeof(subset_t) * (2*num_vertices -1))
     318    tab[0] = 0x7ffffffful >> (31 - num_vertices)
     319
     320    cdef int beg = 0
     321    cdef int end = 1
     322
     323    while beg != end:
     324
     325        s = tab[beg]
     326        beg += 1
     327
     328        if cslots[s] == 0:
     329            continue
     330
     331        g.add_edge(bitset_to_vertex_set(s), bitset_to_vertex_set(s&~cslots[s]))
     332        g.add_edge(bitset_to_vertex_set(s), bitset_to_vertex_set(cslots[s]))
     333
     334        tab[end] = s&~cslots[s]
     335        end += 1
     336        tab[end] = cslots[s]
     337        end += 1
     338
     339    sage_free(tab)
     340    return g
     341
     342cdef bitset_to_vertex_set(subset_t s):
     343    """
     344    Returns as a Set object the set corresponding to the given subset_t
     345    variable.
     346    """
     347    from sage.rings.integer import Integer
     348    from sage.sets.set import Set
     349    from sage.misc.bitset import FrozenBitset
     350    return Set(list(FrozenBitset((Integer(<unsigned int> s).binary())[::-1])))
  • new file sage/graphs/graph_decompositions/rankwidth/README

    diff --git a/sage/graphs/graph_decompositions/rankwidth/README b/sage/graphs/graph_decompositions/rankwidth/README
    new file mode 100644
    - +  
     1This 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.
     2On today's computers it works quite well up to graph sizes of about 28 nodes.
     3
     4Assuming that libigraph is installed it can be compiled using e.g. "gcc -O2 --std=c99 -pedantic rw.c simplerw.c -ligraph".
     5
     6It 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".
     7
     8rw.h / rw.c provide a general mechanism for calculating rank-width and rank-decompositions. simplerw.c provides a simple command-line interface.
     9
  • new file sage/graphs/graph_decompositions/rankwidth/__init__.py

    diff --git a/sage/graphs/graph_decompositions/rankwidth/__init__.py b/sage/graphs/graph_decompositions/rankwidth/__init__.py
    new file mode 100644
    - +  
     1# This file is not empty !
  • new file sage/graphs/graph_decompositions/rankwidth/rw.c

    diff --git a/sage/graphs/graph_decompositions/rankwidth/rw.c b/sage/graphs/graph_decompositions/rankwidth/rw.c
    new file mode 100644
    - +  
     1// Calculate rank-width and rank-decompositions of graphs.
     2
     3// Philipp Klaus Krause, philipp@informatik.uni-frankfurt.de, pkk@spth.de, 2009 - 2012
     4// Copyright (c) 2009-2012 Philipp Klaus Krause
     5
     6// This program is free software; you can redistribute it and/or modify
     7// it under the terms of the GNU General Public License as published by
     8// the Free Software Foundation; either version 2 of the License, or
     9// (at your option) any later version.
     10//
     11// This program is distributed in the hope that it will be useful,
     12// but WITHOUT ANY WARRANTY; without even the implied warranty of
     13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     14// GNU General Public License for more details.
     15//
     16// You should have received a copy of the GNU General Public License along
     17// with this program; if not, write to the Free Software Foundation, Inc.,
     18// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
     19   
     20#include <stdlib.h>
     21
     22#include "rw.h"
     23
     24static uint_fast8_t subset_size;
     25
     26static uint_fast8_t num_vertices;
     27
     28subset_t *adjacency_matrix;
     29
     30// This array is meant to contain max{f(Y), w(f, Y)} at the position corresponding to Y.
     31static uint_fast8_t *slots;
     32subset_t *cslots;
     33
     34uint_fast8_t cut_rank(const subset_t s)
     35{
     36        subset_t am[num_vertices];
     37        subset_t x, y;
     38        uint_fast8_t i, j;
     39        uint_fast8_t rank = 0;
     40
     41        for(i = 0; i < num_vertices; i++)
     42                am[i] =  (s & (1ul << i)) ? 0 : (adjacency_matrix[i] & s);
     43
     44        for(i = 0; i < num_vertices; i++)
     45        {
     46                y = 0;
     47                for(j = rank; j < num_vertices; j++)
     48                {
     49                        x = am[j];
     50                        if(x & 1)
     51                        {
     52                                if(!y)
     53                                {
     54                                        y = x;
     55                                        am[j] = am[rank++];
     56                                        continue;
     57                                }
     58                                x ^= y;
     59                        }
     60                        am[j] = x >> 1;
     61                }
     62        }
     63
     64        return(rank);
     65}
     66
     67// Compute the binomial coefficient C(n, k).
     68// Time complexity: O(n) * complexity of integer division.
     69// Could be replaced by a lookup table of size O(n^2), but since this is not the bopttleneck, it doesn't matter.
     70subset_t binomial_coefficient(uint_fast8_t n, uint_fast8_t k)
     71{
     72        uint_fast8_t i, delta, max;
     73        subset_t c;
     74
     75        if(n < k)
     76                return(0);
     77
     78        if(n == k)
     79                return(1);
     80
     81        if(k < n - k)
     82        {
     83                delta = n - k;
     84                max = k;
     85        }
     86        else
     87        {
     88                delta = k;
     89                max = n - k;
     90        }
     91
     92        c = delta + 1;
     93
     94        for(i = 2; i <= max; i++)
     95                c = (c * (delta + i)) / i;
     96
     97        return(c);
     98}
     99
     100// Convert unsigned integer from combinadic to machine representation.
     101subset_t comb_to_int(subset_t c)
     102{
     103        uint_fast8_t k, n;
     104        subset_t i;
     105
     106        i = 0;
     107        for(k = 0, n = 1; k < num_vertices; k++, c >>= 1)
     108                if(c & 1)
     109                        i += binomial_coefficient(k, n++);
     110        return(i);
     111}
     112
     113// Return largest value v where v < a and  Choose(v,b) <= x.
     114static uint_fast8_t largest_v(uint_fast8_t a, uint_fast8_t b, subset_t x)
     115{
     116        uint_fast8_t v = a - 1;
     117           
     118        while(binomial_coefficient(v,b) > x)
     119                v--;
     120
     121        return(v);
     122}
     123
     124// Convert unsigned integer from machine representation to combinadic.
     125static subset_t int_to_comb(subset_t i)
     126{
     127        uint_fast8_t j, v;
     128        uint_fast8_t a = num_vertices;
     129        uint_fast8_t b = subset_size;
     130        subset_t c = 0;
     131
     132        for(j = 0; j < subset_size; j++, b--)
     133        {
     134                v = largest_v(a, b, i);
     135                i = i - binomial_coefficient(v, b);
     136                a = v;
     137                c |= (1ul << v);
     138        }
     139
     140        return(c);
     141}
     142
     143// Masked increment.
     144static subset_t subset_inc(subset_t v, subset_t mask)
     145{
     146        return((v - mask) & mask);
     147}
     148
     149// 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.
     150// This is where most of the time is spent.
     151uint_fast8_t width(subset_t s)
     152{
     153        uint_fast8_t w = UINT_FAST8_MAX, v, v1, v2;
     154        subset_t ss;
     155        subset_t cs;
     156        for(ss = subset_inc(0, s); ss != s; ss = subset_inc(ss, s))
     157        {
     158                v1 = slots[ss];
     159                v2 = slots[s & ~ss];
     160                v = v1 > v2 ? v1 : v2;
     161                if(v < w)
     162                {
     163                        w = v;
     164                        cs = ss;
     165                }
     166        }
     167        cslots[s] = cs;
     168        return(w);
     169}
     170
     171void fill_slot(subset_t i)
     172{
     173        uint_fast8_t v, w;
     174        subset_t s = int_to_comb(i);
     175        v = cut_rank(s);
     176        w = width(s);
     177        slots[s] = v > w ? v : w;
     178}
     179
     180void calculate_level(uint_fast8_t ss)
     181{
     182        uint_fast8_t i;
     183
     184        subset_size = ss;
     185
     186        if(subset_size == 0)
     187                slots[0] = 0;
     188        else if(subset_size == 1)
     189                for(i = 0; i < num_vertices; i++)
     190                {
     191                        slots[1ul << i] = cut_rank(1ul << i);
     192                        cslots[1ul << i] = 0x0;
     193                }
     194        else
     195        {
     196                subset_t i;
     197                const subset_t end = binomial_coefficient(num_vertices, subset_size);
     198                for(i = 0; i < end; i++)
     199                        fill_slot(i);
     200        }
     201}
     202
     203void calculate_all(void)
     204{
     205        uint_fast8_t i;
     206
     207        for(i = 0; i < num_vertices; i++)
     208        {
     209                slots[1ul << i] = cut_rank(1ul << i);
     210                cslots[1ul << i] = 0x0;
     211        }
     212
     213        for(subset_size = 2; subset_size <= num_vertices; subset_size++)
     214        {
     215                subset_t i;
     216                const subset_t end = binomial_coefficient(num_vertices, subset_size);
     217                for(i = 0; i < end; i++)
     218                        fill_slot(i);
     219        }
     220}
     221
     222int init_rw(uint_fast8_t n)
     223{
     224        // If sizeof(uint_fast8_t) * (1ul << n) overflows, it wraps around to 0, since size_t and unsigned long are unsigned integer types.
     225        if(n > MAX_VERTICES || n && !(sizeof(uint_fast8_t) * (1ul << n)))
     226                return(-1);
     227
     228        num_vertices = n;
     229        adjacency_matrix = malloc(sizeof(subset_t) * n);
     230        slots = malloc(sizeof(uint_fast8_t) * (1ul << n));
     231        cslots = 0;
     232        return((adjacency_matrix && slots) ? 0 : -1);
     233}
     234
     235int init_rw_dec(uint_fast8_t n)
     236{
     237        // If sizeof(subset_t) * (1ul << n) overflows, it wraps around to 0, since size_t and unsigned long are unsigned integer types.
     238        if(n && !(sizeof(subset_t) * (1ul << n)))
     239                return(-1);
     240
     241        if(init_rw(n))
     242                return(-1);
     243        cslots = malloc(sizeof(subset_t) * (1ul << n));
     244        return(cslots ? 0 : -1);
     245}
     246
     247void destroy_rw(void)
     248{
     249        free(cslots);
     250        free(slots);
     251        free(adjacency_matrix);
     252}
     253
     254uint_fast8_t get_rw(void)
     255{
     256        return(slots[0xfffffffful >> (32 - num_vertices)]);
     257}
     258
  • new file sage/graphs/graph_decompositions/rankwidth/rw.h

    diff --git a/sage/graphs/graph_decompositions/rankwidth/rw.h b/sage/graphs/graph_decompositions/rankwidth/rw.h
    new file mode 100644
    - +  
     1// Calculate rank-width and rank-decompositions of graphs.
     2
     3// Philipp Klaus Krause, philipp@informatik.uni-frankfurt.de, pkk@spth.de, 2009 - 2011
     4// Copyright (c) 2009-2011 Philipp Klaus Krause
     5
     6// This program is free software; you can redistribute it and/or modify
     7// it under the terms of the GNU General Public License as published by
     8// the Free Software Foundation; either version 2 of the License, or
     9// (at your option) any later version.
     10//
     11// This program is distributed in the hope that it will be useful,
     12// but WITHOUT ANY WARRANTY; without even the implied warranty of
     13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     14// GNU General Public License for more details.
     15//
     16// You should have received a copy of the GNU General Public License along
     17// with this program; if not, write to the Free Software Foundation, Inc.,
     18// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
     19
     20// 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.
     21// On 2009's computers it works quite well up to graph sizes of about 28 nodes.
     22// 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.
     23
     24#include <stdint.h>
     25
     26// 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).
     27// uint_leastN_t is faster than uint_fastN_t here, since the bottleneck is cache misses.
     28typedef uint_least32_t subset_t;
     29#define MAX_VERTICES 32
     30
     31// Input graph.
     32//extern subset_t adjacency_matrix[NUM_VERTICES];
     33extern subset_t *adjacency_matrix;
     34
     35// Output rank-decomposition
     36extern subset_t *cslots;
     37
     38// Initialization (for getting rank-width only). Returns 0 on success.
     39//int init_rw(uint_fast8_t n);
     40
     41// Initialization (for getting both rank-width and rank-decomposition). Returns 0 on success.
     42int init_rw_dec(uint_fast8_t n);
     43
     44// Free memory allocated during initialization.
     45void destroy_rw(void);
     46
     47// Calculate everything. May take some time.
     48void calculate_all(void);
     49
     50// Calculate a single level only
     51void calculate_level(uint_fast8_t subset_size);
     52
     53// Get the rank-width.
     54uint_fast8_t get_rw(void);
     55
  • new file sage/graphs/graph_decompositions/rankwidth/simplerw.c

    diff --git a/sage/graphs/graph_decompositions/rankwidth/simplerw.c b/sage/graphs/graph_decompositions/rankwidth/simplerw.c
    new file mode 100644
    - +  
     1// A simple program for calculating rank-width and rank-decompositions.
     2// Compile using gcc -O2 --std=c99 -pedantic rw.c simplerw.c -ligraph
     3
     4// Philipp Klaus Krause, philipp@informatik.uni-frankfurt.de, pkk@spth.de, 2009 - 2011
     5// Copyright (c) 2009-2011 Philipp Klaus Krause
     6
     7// This program is free software; you can redistribute it and/or modify
     8// it under the terms of the GNU General Public License as published by
     9// the Free Software Foundation; either version 2 of the License, or
     10// (at your option) any later version.
     11//
     12// This program is distributed in the hope that it will be useful,
     13// but WITHOUT ANY WARRANTY; without even the implied warranty of
     14// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     15// GNU General Public License for more details.
     16//
     17// You should have received a copy of the GNU General Public License along
     18// with this program; if not, write to the Free Software Foundation, Inc.,
     19// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
     20
     21#include <stdio.h>
     22#include <stdlib.h>
     23#include <string.h>
     24
     25#include <igraph/igraph.h>
     26
     27#include "rw.h"
     28
     29static int num_vertices = 18;
     30
     31static uint_fast32_t bitmask(int i)
     32{
     33        return(1ul << i);
     34}
     35
     36static void set_am(int i, int j, int val)
     37{
     38        adjacency_matrix[i] &= ~bitmask(j);
     39        adjacency_matrix[j] &= ~bitmask(i);
     40        if(val)
     41        {
     42                adjacency_matrix[i] |= bitmask(j);
     43                adjacency_matrix[j] |= bitmask(i);
     44        }
     45}
     46
     47static void tab(int l)
     48{
     49        while(l--)
     50                printf("\t");
     51}
     52
     53void print_rank_dec(subset_t s, int l)
     54{
     55        tab(l);
     56        printf("cslot: %x\n", (unsigned)s);
     57        if(cslots[s] == 0)
     58                return;
     59        print_rank_dec(cslots[s], l + 1);
     60        print_rank_dec(s & ~cslots[s], l + 1);
     61}
     62
     63/*int random_graph(void)
     64{
     65        int i, j;
     66
     67        if(init_rw_dec(num_vertices))
     68                return(-1);
     69
     70        srand(23);
     71        for(i = 0; i < num_vertices; i++)
     72                for(j = 0; j < num_vertices; j++)
     73                        set_am(i, j, rand() % 2);
     74
     75        for(i = 0; i < num_vertices; i++)
     76                printf("Adj. mat.: %x\n", (unsigned)(adjacency_matrix[i]));
     77
     78        return(0);
     79}*/
     80
     81int read_graph(const char *format, const char * filename)
     82{
     83        int ret, i, j;
     84        FILE *file;
     85        igraph_t igraph;
     86        igraph_matrix_t imatrix;
     87        igraph_bool_t isimple;
     88
     89        // Open file
     90        if(!(file = fopen(filename, "r")))
     91        {
     92                printf("Failed to open file %s.\n", filename);
     93                return(-1);
     94        }
     95       
     96        // Read graph from file
     97        if(!strcmp("--edgelist", format))
     98                ret = igraph_read_graph_edgelist(&igraph, file, 0, 0);
     99        else if(!strcmp("--graphdb", format))
     100                ret = igraph_read_graph_graphdb(&igraph, file, 0);
     101        else if(!strcmp("--graphml", format))
     102                ret = igraph_read_graph_graphml(&igraph, file, 0);
     103        else if(!strcmp("--gml", format))
     104                ret = igraph_read_graph_gml(&igraph, file);
     105        else if(!strcmp("--pajek", format))
     106                ret = igraph_read_graph_pajek(&igraph, file);
     107        else
     108        {
     109                printf("Unknown file format %s.\n", format);
     110                fclose(file);
     111                return(-1);
     112        }
     113
     114        fclose(file);
     115
     116        if(ret)
     117        {
     118                printf("Failed to read a graph from file %s.\n", filename);
     119                return(-1);
     120        }
     121
     122        // Check for undirectedness
     123        if(igraph_is_directed(&igraph) || igraph_is_simple(&igraph, &isimple) || !isimple)
     124        {
     125                printf("Input is not a simple, undirected graph from file %s.\n", filename);
     126                igraph_destroy(&igraph);
     127                return(-1);
     128        }
     129
     130        // Get adjacency matrix
     131        if(igraph_matrix_init(&imatrix, 1, 1))
     132        {
     133                igraph_destroy(&igraph);
     134                return(-1);
     135        }
     136        igraph_get_adjacency(&igraph, &imatrix, IGRAPH_GET_ADJACENCY_BOTH);
     137        igraph_destroy(&igraph);
     138        if(igraph_matrix_nrow(&imatrix) > MAX_VERTICES)
     139        {
     140                igraph_matrix_destroy(&imatrix);
     141                printf("Graph too large.\n");
     142                return(-1);
     143        }
     144
     145        num_vertices = igraph_matrix_nrow(&imatrix);
     146       
     147        if(init_rw_dec(num_vertices))
     148        {
     149                destroy_rw();
     150                printf("Not enough memory.\n");
     151                return(-1);
     152        }
     153
     154        for(i = 0; i < num_vertices; i++)
     155                for(j = 0; j < num_vertices; j++)
     156                        set_am(i, j, MATRIX(imatrix, i, j));
     157
     158        igraph_matrix_destroy(&imatrix);
     159       
     160        return(ret ? -1 : 0);
     161}
     162
     163int main(int argc, char *argv[])
     164{
     165        int i;
     166
     167        if(argc /*==*/ <= 2 || argc > 4)
     168        {
     169                printf("Usage: rw [--format filename]\n");
     170                printf("Supported formats: edgelist, graphdb, graphml, gml, pajek.\n");
     171                return(-1);
     172        }
     173
     174        if(/*argc <= 1 ? random_graph() :*/ read_graph(argv[1], argv[2]))
     175        {
     176                printf("Failed to create input graph.\n");
     177                return(-1);
     178        }
     179
     180        printf("%d vertices.\n", (int)num_vertices);
     181
     182        for(i = 0; i <= num_vertices; i++)
     183        {
     184                printf("Calculating for subsets of size %d.\n", i);
     185                calculate_level(i);
     186        }
     187
     188        printf("rank-width: %d\n", (int)get_rw());
     189
     190        print_rank_dec(0x7ffffffful >> (31 - num_vertices), 0);
     191
     192        destroy_rw();
     193
     194        return(0);
     195}
     196
  • sage/graphs/modular_decomposition/__init__.py

    diff --git a/sage/graphs/modular_decomposition/__init__.py b/sage/graphs/modular_decomposition/__init__.py
    a b  
    11# This file is not empty !
     2
     3import sage.graphs.modular_decomposition.modular_decomposition
     4
  • sage/graphs/modular_decomposition/modular_decomposition.pyx

    diff --git a/sage/graphs/modular_decomposition/modular_decomposition.pyx b/sage/graphs/modular_decomposition/modular_decomposition.pyx
    a b  
     1r"""
     2Modular decomposition
     3"""
     4
    15#####################################################
    26# The following code is mainly a Cythonized
    37# copy of code found in src/random.c and src/dm.c
     
    2630
    2731cpdef modular_decomposition(g):
    2832    r"""
    29     Returns a modular decomposition of the given graph
     33    Returns a modular decomposition of the given graph.
    3034
    3135    INPUT:
    3236
     
    4246            * ``"Prime"``
    4347            * ``"Serie"``
    4448
    45         * The list of submodules (as list of pairs ``(type, list)``, recursively...)
    46           or the vertex's name if the module is a singleton.         
     49        * The list of submodules (as list of pairs ``(type, list)``,
     50          recursively...)  or the vertex's name if the module is a
     51          singleton.
    4752
    4853    .. NOTE::
    4954
    50     As this fuction could be used by efficient C routines, the vertices returned are
    51     not labels but identifiants from ``[0, ..., g.order()-1]``
    52 
    53 
     55        As this fuction could be used by efficient C routines, the
     56        vertices returned are not labels but identifiants from ``[0,
     57        ..., g.order()-1]``
    5458
    5559    ALGORITHM:
    5660
    5761    This function uses a C implementation of a 2-step algorithm
    58     implemented by Fabien de Montgolfier [FMDec]_ :
     62    implemented by Fabien de Montgolfier [FMDecb]_ :
    5963
    60         * Computation of a factorizing permutation [HabibViennot1999]_.
     64        * Computation of a factorizing permutation [HabibViennot1999b]_.
    6165
    62         * Computation of the tree itself [CapHabMont02]_.
    63 
    64     REFERENCE:
    65 
    66     .. [FMDec] Fabien de Montgolfier
    67       http://www.liafa.jussieu.fr/~fm/algos/index.html
    68 
    69     .. [HabibViennot1999] Michel Habib, Christiphe Paul, Laurent Viennot
    70       Partition refinement techniques: An interesting algorithmic tool kit
    71       International Journal of Foundations of Computer Science
    72       vol. 10 n2 pp.147--170, 1999
    73 
    74     .. [CapHabMont02] C. Capelle, M. Habib et F. de Montgolfier
    75       Graph decomposition and Factorising Permutations
    76       Discrete Mathematics and Theoretical Computer Sciences, vol 5 no. 1 , 2002.
     66        * Computation of the tree itself [CapHabMont02b]_.
    7767
    7868    EXAMPLES:
    7969
     
    9686        sage: g.add_edge(0,6)
    9787        sage: modular_decomposition(g)
    9888        ('Serie', [0, ('Parallel', [5, ('Serie', [1, 4, 3, 2]), 6])])
     89
     90
     91    REFERENCES:
     92
     93    .. [FMDecb] Fabien de Montgolfier
     94      http://www.liafa.jussieu.fr/~fm/algos/index.html
     95
     96    .. [HabibViennot1999b] Michel Habib, Christiphe Paul, Laurent Viennot
     97      Partition refinement techniques: An interesting algorithmic tool kit
     98      International Journal of Foundations of Computer Science
     99      vol. 10 n2 pp.147--170, 1999
     100
     101    .. [CapHabMont02b] C. Capelle, M. Habib et F. de Montgolfier
     102      Graph decomposition and Factorising Permutations
     103      Discrete Mathematics and Theoretical Computer Sciences, vol 5 no. 1 , 2002.
    99104    """
    100105    cdef c_graphe G
    101106    cdef c_adj * a
  • setup.py

    diff --git a/setup.py b/setup.py
    a b  
    905905                     'sage.graphs.base',
    906906                     'sage.graphs.modular_decomposition',
    907907                     'sage.graphs.graph_decompositions',
     908                     'sage.graphs.graph_decompositions',
    908909                     
    909910                     'sage.groups',
    910911                     'sage.groups.abelian_gps',