r"""
Sage Sandpiles

Functions and classes for mathematical sandpiles.

Version: 1.51

CHANGES
-- fixed bug with calculating linear systems and effective divisors

AUTHOR:
    -- David Perkinson (2008-12-27): initial version

    -- David Perkinson (2009): many undocumented improvements
    
    -- David Perkinson (2009-07-15): switched to using config_to_list instead
       of .values(), thus fixing a few bugs when not using integer labels for
       vertices.

EXAMPLES:

A weighted directed graph given as a Python dictionary:

    sage: g = {0: {},
	       1: {0: 1, 2: 1, 3: 1},
	       2: {1: 1, 3: 1, 4: 1},
	       3: {1: 1, 2: 1, 4: 1},
	       4: {2: 1, 3: 1}}

The associated sandpile with 0 chosen as the sink::

    sage: S = Sandpile(g,0)

A picture of the graph::
   
    sage: S.show()

The relevant Laplacian matrices::

    sage: S.laplacian()
    [ 0  0  0  0  0]
    [-1  3 -1 -1  0]
    [ 0 -1  3 -1 -1]
    [ 0 -1 -1  3 -1]
    [ 0  0 -1 -1  2]
    sage: S.reduced_laplacian()
    [ 3 -1 -1  0]
    [-1  3 -1 -1]
    [-1 -1  3 -1]
    [ 0 -1 -1  2]

The number of elements of the sandpile group for S::

    sage: S.group_order()
    8

The structure of the sandpile group::

    sage: S.elementary_divisors()
    [1, 1, 1, 8]

The elements of the sandpile group for S::
          
    sage: S.recurrents()
    [{1: 2, 2: 2, 3: 2, 4: 1},
     {1: 2, 2: 2, 3: 2, 4: 0},
     {1: 2, 2: 1, 3: 2, 4: 0},
     {1: 2, 2: 2, 3: 0, 4: 1},
     {1: 2, 2: 0, 3: 2, 4: 1},
     {1: 2, 2: 2, 3: 1, 4: 0},
     {1: 2, 2: 1, 3: 2, 4: 1},
     {1: 2, 2: 2, 3: 1, 4: 1}]

The maximal stable element (2 grains of sand on vertices 1, 2, and 3, and 1
grain of sand on vertex 4::

    sage: S.max_stable()
    {1: 2, 2: 2, 3: 2, 4: 1}
    sage: S.config_to_list(S.max_stable())
    [2, 2, 2, 1]

The identity of the sandpile group for S::

    sage: S.identity()
    {1: 2, 2: 2, 3: 2, 4: 0}

Some group operations::

    sage: c = S.add(S.max_stable(),S.identity())
    sage: c
    {1: 4, 2: 4, 3: 4, 4: 1}

Firing an unstable vertex returns resulting configuration::

    sage: S.fire_vertex(1,c)
    {1: 1, 2: 5, 3: 5, 4: 1}
    sage: c
    {1: 1, 2: 5, 3: 5, 4: 1}

Fire all currently unstable vertices, returning the resulting configuration and
the firing vector::

    sage: S.stabilize_one_step(c)
    [{1: 3, 2: 3, 3: 3, 4: 3}, {1: 0, 2: 1, 3: 1, 4: 0}]
    sage: c
    {1: 3, 2: 3, 3: 3, 4: 3}

Finish stabilizing c, returning the resulting configuration and the firing
vector::

    sage: S.stabilize(c)
    [{1: 2, 2: 2, 3: 2, 4: 1}, {1: 5, 2: 7, 3: 7, 4: 8}]
    sage: c
    {1: 2, 2: 2, 3: 2, 4: 1}

Similarly,

    sage: S.stably_add(S.max_stable(),S.identity())
    [{1: 2, 2: 2, 3: 2, 4: 1}, {1: 6, 2: 8, 3: 8, 4: 8}]

The number of superstable configurations of each degree::

    sage: S.hilbert_function()
    [1, 4, 8]
    sage: S.postulation()
    2

the saturated, homogeneous sandpile ideal

    sage: S.ideal()    

    x_1-x_0,
    x_3*x_2-x_0^2,
    x_4^2-x_0^2,
    x_2^3-x_4*x_3*x_0,
    x_4*x_2^2-x_3^2*x_0,
    x_3^3-x_4*x_2*x_0,
    x_4*x_3^2-x_2^2*x_0

its minimal free resolution

    sage: S.resolution()
    'R <-- R^7 <-- R^15 <-- R^13 <-- R^4'

and its Betti numbers::

    sage: S.betti()
               0     1     2     3     4
    ------------------------------------
        0:     1     1     -     -     -
        1:     -     2     2     -     -
        2:     -     4    13    13     4
    ------------------------------------
    total:     1     7    15    13     4

Distribution of avalanche sizes::

    sage: S = Sandpile(grid(10,10),'sink')
    sage: S.set_config(S.max_stable())
    sage: a = [sum(S.stabilize(S.add_random())[1].values()) for i in range(1000)]
    sage: p = list_plot([[log(i+1),log(a.count(i))] for i in [0..max(a)] if a.count(i)])
    sage: t = text("Distribution of avalanche sizes", (2,2), rgbcolor=(1,0,0))
    sage: show(p+t)
"""

# set the following if 4ti2 is installed
path_to_zsolve = '/home/davidp/math/sandpile/sage/sage-sandpile1.5/4ti2/linux_x86/'

#*****************************************************************************
#       Copyright (C) 2009 David Perkinson <davidp@reed.edu>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#                  http://www.gnu.org/licenses/
#*****************************************************************************

class Sandpile(DiGraph):
    """
    Class for Dhar's abelian sandpile model.
    """

    def __init__(self, g, sink):
        r"""
        Create a sandpile.

        INPUT: 
	
	 - ``g`` - dict for directed multgraph (see NOTES) edges weighted by
	   nonnegative integers

	 - ``sink`` - A sink vertex.  Any outgoing edges from the designated
	   sink are ignored for the purposes of stabilization.  It is assumed
	   that every vertex has a directed path into the sink.

        OUTPUT:

	sandpile

	EXAMPLES:

        Here ``g`` represents a square with directed, multiple edges with three
        vertices, ``a``, ``b``, ``c``, and ``d``.  The vertex ``a`` has
        outgoing edges to itself (weight 2), to vertex ``b`` (weight 1), and
        vertex ``c`` (weight 3), for example.

        ::

            sage: g = {'a': {'a':2, 'b':1, 'c':3}, 'b': {'a':1, 'd':1},
                       'c': {'a':1,'d': 1}, 'd': {'b':1, 'c':1}}
            sage: G = Sandpile(g,'d')

	Here is a square with unweighted edges.  In this example, the graph is
	also undirected.

	::

            sage: g = {0:[1,2], 1:[0,3], 2:[0,3], 3:[1,2]}
            sage: G = Sandpile(g,3)
       
        NOTES::

	Loops are allowed.  There are two admissible formats, both of which are
	dictionaries whose keys are the vertex names.  In one format, the
	values are dictionaries with keys the names of vertices which are the
	tails of outgoing edges and with values the weights of the edges.  In
	the other format, the values are lists of names of vertices which are
	the tails of the outgoing edges.  All weights are then automatically
	assigned the value 1.
        """
        # preprocess a graph, if necessary
        if type(g) == dict and type(g.values()[0]) == dict:
            pass # this is the default format
        elif type(g) == dict and type(g.values()[0]) == list:
            processed_g = {}
            for k in g.keys():
                temp = {}
                for vertex in g[k]:
                    temp[vertex] = 1
                processed_g[k] = temp
            g = processed_g
        elif type(g) == sage.graphs.graph.Graph:
            processed_g = {}
            for v in g.vertices():
                edges = {}
                for n in g.neighbors(v):
                    if type(g.edge_label(v,n)) == type(1) and g.edge_label(v,n) >=0:
                        edges[n] = g.edge_label(v,n)
                    else:
                        edges[n] = 1
                processed_g[v] = edges
            g = processed_g
        elif type(g) == sage.graphs.graph.DiGraph:
            processed_g = {}
            for v in g.vertices():
                edges = {}
                for n in g.successors(v):
                    if (type(g.edge_label(v,n)) == type(1)
                        and g.edge_label(v,n)>=0):
                        edges[n] = g.edge_label(v,n)
                    else:
                        edges[n] = 1
                processed_g[v] = edges
            g = processed_g
        else:
            raise SyntaxError, g

        # create digraph and initialize some variables
        DiGraph.__init__(self,g,weighted=True)
        self._dict = deepcopy(g)
        self._sink = sink  # key for sink
        self._sink_ind = self.vertices().index(sink)
        self._nonsink_vertices = deepcopy(self.vertices())
        del self._nonsink_vertices[self._sink_ind]
        # compute laplacians:
        self._laplacian = self.laplacian_matrix()
        temp = range(self.num_verts())
        del temp[self._sink_ind]
        self._reduced_laplacian = self._laplacian[temp,temp]
        # the current configuration is initially the all-zero configuration
        self._config = {}
        self.reset_config(self._config)

    def __getattr__(self, name):
        """
        Set certain variables only when called.
        """
        if not self.__dict__.has_key(name):
            if name == '_max_stable':
                self._set_max_stable()
                return deepcopy(self.__dict__[name])
            elif name == '_out_degrees':
                self._set_out_degrees()
                return deepcopy(self.__dict__[name])
            elif name == '_in_degrees':
                self._set_in_degrees()
                return deepcopy(self.__dict__[name])
            elif name == '_burning_config' or name == '_burning_script':
                self._set_burning_config()
                return deepcopy(self.__dict__[name])
            elif name == '_identity':
                self._set_identity()
                return deepcopy(self.__dict__[name])
            elif name == '_recurrents':
                self._set_recurrents()
                return deepcopy(self.__dict__[name])
            elif name == '_superstables':
                self._set_superstables()
                return deepcopy(self.__dict__[name])
            elif name == '_group_order':
                self.__dict__[name] = det(self._reduced_laplacian.dense_matrix())
                return self.__dict__[name]
            elif name == '_elementary_divisors':
                self._set_elementary_divisors()
                return deepcopy(self.__dict__[name])
            elif (name == '_postulation' or name == '_first_diffs_hilb' 
                   or name == '_hilbert_function'):
                self._set_hilbert_function()
                return deepcopy(self.__dict__[name])
            elif (name == '_ring' or name == '_variables' 
                   or name == '_unsaturated_ideal'):
                self._set_ring()
                return self.__dict__[name]
            elif name == '_ideal':
                self._set_ideal()
                return self.__dict__[name]
            elif name == '_resolution':
                self._set_resolution()
                return self.__dict__[name]
            elif name == '_groebner':
                self._set_groebner()
                return self.__dict__[name]
            elif name == '_betti':
                self._set_betti()
                return self.__dict__[name]
            elif name == '_points':
                self._set_points()
                return self.__dict__[name]
            else:
                raise AttributeError, name

    def version(self):
        r"""
        Returns the version number of Sage Sandpiles.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	string


	EXAMPLES::

            sage: S = sandlib('generic')
            sage: S.version()
        """
        print 'Sage Sandpiles Version 1.51'

    def dict(self):
        r"""
        Returns a dictionary of dictionaries representing a directed graph.

        INPUT: 
	
	None

        OUTPUT: 
	
	dict


	EXAMPLES::

            sage: G = sandlib('generic')
            sage: G.dict()
	    {0: {},
	     1: {0: 1, 3: 1, 4: 1},
	     2: {0: 1, 3: 1, 5: 1},
	     3: {2: 1, 5: 1},
	     4: {1: 1, 3: 1},
	     5: {2: 1, 3: 1}}
	    sage: G.sink()
	    0
        """
        return deepcopy(self._dict)

    def sink(self):
        r"""
        Returns the identifier for the sink vertex.
	
        INPUT:
	
	None
	
        OUTPUT: 
	Object (name for the sink vertex)

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.sink()
	    0
	    sage: H = Sandpile(grid(2,2),'sink')
	    sage: H.sink()
	    'sink'
	    sage: type(H.sink())
	    <type 'str'>
        """
        return self._sink

    def laplacian(self):
        r"""
        Returns the Laplacian matrix of the graph.

        INPUT: 
	
	None

        OUTPUT: 
	
	matrix


	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.laplacian()
	    [ 0  0  0  0  0  0]
	    [-1  3  0 -1 -1  0]
	    [-1  0  3 -1  0 -1]
	    [ 0  0 -1  2  0 -1]
	    [ 0 -1  0 -1  2  0]
	    [ 0  0 -1 -1  0  2]
        """
        return deepcopy(self._laplacian)

    def reduced_laplacian(self):
        r"""
	Returns the reduced Laplacian matrix of the graph. 
	
        INPUT: 
	
	None
        
	OUTPUT: 
	
	matrix


	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.laplacian()
	    [ 0  0  0  0  0  0]
	    [-1  3  0 -1 -1  0]
	    [-1  0  3 -1  0 -1]
	    [ 0  0 -1  2  0 -1]
	    [ 0 -1  0 -1  2  0]
	    [ 0  0 -1 -1  0  2]
	    sage: G.reduced_laplacian()
	    [ 3  0 -1 -1  0]
	    [ 0  3 -1  0 -1]
	    [ 0 -1  2  0 -1]
	    [-1  0 -1  2  0]
	    [ 0 -1 -1  0  2]

	NOTES:

	This is the Laplacian matrix with the row and column indexed by the
	sink vertex removed.
        """
        return deepcopy(self._reduced_laplacian)

    def group_order(self):
        r"""
	Returns the size of the sandpile group. 

        INPUT: 
	
	None
        
	OUTPUT: 
	
	int

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.group_order()
	    15
        """
        return self._group_order

    def element_order(self, config=None):
        r"""
        Returns the order of the recurrent element equivalent to ``config``.

        INPUT:

        config - configuration

        OUTPUT:

        integer

        EXAMPLES::

            sage: S = sandlib('generic')
            sage: [S.element_order(r) for r in S.recurrents()]
            [1, 3, 15, 5, 5, 15, 5, 15, 3, 15, 15, 5, 15, 15, 15]

        """
        if config is None:
            config = self._config
        v = vector(self.config_to_list(config))
        w = v*self.reduced_laplacian().dense_matrix()^(-1)
        return lcm([denominator(i) for i in w])

    def _set_max_stable(self):
        r"""
        Initialize the variable holding the maximal stable configuration.

	INPUT:

	None

	OUTPUT:

        NONE

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G._set_max_stable()
	"""
        self._max_stable = {}
        for v in self._nonsink_vertices:
            self._max_stable[v] = self.out_degree(v)-1

    def max_stable(self):
        r"""
        Returns the maximal stable configuration.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	dict (the maximal stable configuration)


	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.max_stable()
	    {1: 2, 2: 2, 3: 1, 4: 1, 5: 1}
        """
        return deepcopy(self._max_stable)

    def _set_out_degrees(self):
        r"""
        Initialize the variable holding the out-degrees.

	INPUT:

	None

	OUTPUT:

        NONE

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G._set_out_degrees()
	"""
        self._out_degrees = {}
        for v in self.vertices():
            self._out_degrees[v] = 0
            for e in self.edges_incident(v):
                self._out_degrees[v] += e[2]

    def out_degree(self, v=None):
        r"""
        Return the out-degree of a vertex or a list of all out-degrees.

	INPUT:

	``v`` (optional) - vertex name

	OUTPUT:

	integer or list of integers

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.out_degree(2)
	    3
	    sage: G.out_degree()
	    {0: 0, 1: 3, 2: 3, 3: 2, 4: 2, 5: 2}
	"""
        if not v is None:
            return self._out_degrees[v]
        else:
            return self._out_degrees

    def _set_in_degrees(self):
        """
        Initialize the variable holding the in-degrees.

	INPUT:

	None

	OUTPUT:

        NONE

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G._set_in_degrees()
        """
        self._in_degrees = {}
        self.reset_config(self._in_degrees)
        self._in_degrees[self._sink] = 0
        for e in self.edges():
            self._in_degrees[e[1]] += e[2]

    def in_degree(self, v=None):
        r"""

        Return the in-degree of a vertex or a list of all in-degrees.

	INPUT:

	``v`` - vertex name or None

	OUTPUT:

	integer or list of integers

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.in_degree(2)
	    2
	    sage: G.in_degree()
	    {0: 2, 1: 1, 2: 2, 3: 4, 4: 1, 5: 2}
	"""
        if not v is None:
            return self._in_degrees[v]
        else:
            return self._in_degrees

    def _set_burning_config(self):
        r"""
        Calculate the minimal burning configuration and its corresponding
	script.

	EXAMPLES::

	    sage: g = {0:{},1:{0:1,3:1,4:1},2:{0:1,3:1,5:1},
		       3:{2:1,5:1},4:{1:1,3:1},5:{2:1,3:1}}
	    sage: G = Sandpile(g,0)
	    sage: G.burning_config()
	    {1: 2, 2: 0, 3: 1, 4: 1, 5: 0}
	    sage: G.config_to_list(G.burning_config())
	    [2, 0, 1, 1, 0]
	    sage: G.burning_script()
	    {1: 1, 2: 3, 3: 5, 4: 1, 5: 4}
	    sage: script = G.config_to_list(G.burning_script())
            sage: script
	    [1, 3, 5, 1, 4]
	    sage: matrix(script)*G.reduced_laplacian()
	    [2 0 1 1 0]
        """
	# TODO: Cythonize!
        d = self._reduced_laplacian.nrows()
        burn = sum(self._reduced_laplacian)
        script=[1]*d  # d 1s
        done = False
        while not done:
            bad = -1
            for i in range(d):
                if burn[i] < 0:
                    bad = i
                    break
            if bad == -1:
                done = True
            else:
                burn += self._reduced_laplacian[bad]
                script[bad]+=1
        b = iter(burn)
        s = iter(script)
        self._burning_config = {}
        self._burning_script = {}
        for v in self._nonsink_vertices:
            self._burning_config[v] = b.next()
            self._burning_script[v] = s.next()

    def burning_config(self):
        r"""
        A minimal burning configuration.

	INPUT:

	None

	OUTPUT:

	dict (configuration)

	EXAMPLES::

	    sage: g = {0:{},1:{0:1,3:1,4:1},2:{0:1,3:1,5:1},
		       3:{2:1,5:1},4:{1:1,3:1},5:{2:1,3:1}}
	    sage: G = Sandpile(g,0)
	    sage: G.burning_config()
	    {1: 2, 2: 0, 3: 1, 4: 1, 5: 0}
	    sage: G.config_to_list(G.burning_config())
	    [2, 0, 1, 1, 0]
	    sage: G.burning_script()
	    {1: 1, 2: 3, 3: 5, 4: 1, 5: 4}
	    sage: script = G.config_to_list(G.burning_script())
            sage: script
	    [1, 3, 5, 1, 4]
	    sage: matrix(script)*G.reduced_laplacian()
	    [2 0 1 1 0]

	NOTES:

	The burning configuration and script are computed using a modified
	version of Speer's script algorithm.  This is a generalization to
	directed multigraphs of Dhar's burning algorithm.

	A *burning configuration* is a nonnegative integer-linear
	combination of the rows of the reduced Laplacian matrix having
	nonnegative entries and such that every vertex has a path from some
	vertex in its support.  The corresponding *burning script* gives
	the integer-linear combination needed to obtain the burning
	configuration.  So if `b` is the burning configuration, `\sigma` is its
	script, and `\tilde{L}` is the reduced Laplacian, then `\sigma\cdot
	\tilde{L} = b`.  The *minimal burning configuration* is the one
	with the minimal script (its components are no larger than the
	components of any other script
	for a burning configuration).

	The following are equivalent for a configuration `c` with burning
	configuration `b` having script `\sigma`:

	 - `c` is recurrent; 
	 - `c+b` stabilizes to `c`;
	 - the firing vector for the stabilization of `c+b` is `\sigma`.
	"""
        return deepcopy(self._burning_config)

    def burning_script(self):
        r"""
        A script for the minimal burning configuration.

	INPUT:

	None

	OUTPUT:

	dict

	EXAMPLES::

	    sage: g = {0:{},1:{0:1,3:1,4:1},2:{0:1,3:1,5:1},
		       3:{2:1,5:1},4:{1:1,3:1},5:{2:1,3:1}}
	    sage: G = Sandpile(g,0)
	    sage: G.burning_config()
	    {1: 2, 2: 0, 3: 1, 4: 1, 5: 0}
	    sage: G.config_to_list(G.burning_config())
	    [2, 0, 1, 1, 0]
	    sage: G.burning_script()
	    {1: 1, 2: 3, 3: 5, 4: 1, 5: 4}
	    sage: script = G.config_to_list(G.burning_script())
            sage: script
	    [1, 3, 5, 1, 4]
	    sage: matrix(script)*G.reduced_laplacian()
	    [2 0 1 1 0]

	NOTES:

	The burning configuration and script are computed using a modified
	version of Speer's script algorithm.  This is a generalization to
	directed multigraphs of Dhar's burning algorithm.

	A *burning configuration* is a nonnegative integer-linear
	combination of the rows of the reduced Laplacian matrix having
	nonnegative entries and such that every vertex has a path from some
	vertex in its support.  The corresponding *burning script* gives the
	integer-linear combination needed to obtain the burning configuration.
	So if `b` is the burning configuration, `s` is its script, and
	`L_{\mathrm{red}}` is the reduced Laplacian, then `s\cdot
	L_{\mathrm{red}}= b`.  The *minimal burning configuration* is the one
	with the minimal script (its components are no larger than the
	components of any other script
	for a burning configuration).

	The following are equivalent for a configuration `c` with burning
	configuration `b` having script `s`:

	 - `c` is recurrent; 
	 - `c+b` stabilizes to `c`;
	 - the firing vector for the stabilization of `c+b` is `s`.
	"""
        return deepcopy(self._burning_script)

    def nonsink_vertices(self):
        r"""
	The names of the nonsink vertices.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.nonsink_vertices()
	    [1, 2, 3, 4, 5]
	"""
        return self._nonsink_vertices

    def unstable(self, config=None):
        r"""
        The list of unstable vertices in ``config``.

	INPUT:
        
        ``config`` (optional) - dict 

	OUTPUT:

        list of vertex names

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: c = G.add(G.max_stable(),G.identity())
	    sage: G.unstable(c)
	    [1, 2, 3, 4]
	    sage: G.stabilize(c);
	    sage: G.unstable(c)
	    []
        """
        if config is None:
            config = self._config
	# the list of unstable vertices in config:
        result = []
        for v in self._nonsink_vertices:
            if config[v] >= self.out_degree(v):
                result.append(v)
        return result

    def config(self):
        r"""
	The current saved configuration.

        INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.set_config(G.max_stable())
	    sage: G.config() == G.max_stable()
	    True
	    sage: G.add(G.max_stable(),G.identity())
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 1}
	    sage: G.config()
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 1}

	NOTES:  
	
	Each sandpile has a saved configuration (initially the zero
	configuration)  The configuration can be set with the ``set_config()``
	method and may be affected by certain other methods, e.g.,
	``stabilize()``, ``add()``, and ``add_random()``---usually those methods
	that output configurations.
	"""
        return deepcopy(self._config)

    def set_config(self, config):
        r"""
        Sets the current configuration to ``config``.

	INPUT:

	``config`` - dict

	OUTPUT:

	None

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.set_config(G.identity())
	    sage: G.config() == G.identity()
	    True
        """
        self._config = deepcopy(config)

    def reset_config(self, config=None):
        r"""
        Turn ``config`` into the zero configuration.

	INPUT:

	``config`` (optional) - dict

	OUTPUT:

	None

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: G.add(G.max_stable(),G.max_stable())
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 2}
	    sage: G.config()
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 2}
	    sage: G.reset_config()
	    sage: G.config()
	    {1: 0, 2: 0, 3: 0, 4: 0, 5: 0}
        """
        if config is None:
            config = self._config
        for v in self._nonsink_vertices:
            config[v] = 0

    def fire_vertex(self, v, config=None):
        r"""
	Fire (topple) a given vertex of a configuration.
	
	INPUT:

	 - ``v`` - vertex name
	 - ``config`` - dict

	OUTPUT:
         
	dict (configuration)

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: c = G.add(G.max_stable(),G.identity())
	    sage: G.unstable(c)
	    [1, 2, 3, 4]
	    sage: c
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 1}
	    sage: G.fire_vertex(1,c)
	    {1: 1, 2: 4, 3: 3, 4: 3, 5: 1}

	NOTES:

	This method fires vertex ``v`` in ``config`` provided ``v`` is a
	nonsink, unstable vertex.  Returns the result (and modifies
	``self.config``).  The vertex ``v`` is fired only once, even if the
	result leaves ``v`` unstable.
        """
        if config is None:
            config = self._config
        if v != self._sink and config[v] >= self.out_degree(v):
            config[v] -= self.out_degree(v)
            for e in self.outgoing_edges(v):
                if e[1] != self._sink:
                    config[e[1]]+=e[2]
            self._config = config
        return config

    def stabilize_one_step(self, config=None):
        r"""
	Fire each unstable vertex in ``config``, returning a list
	``[out_config, firing_vector]`` where ``out_config`` is the modified
	configuration.

	INPUT:

	``config`` (optional) - dict

	OUTPUT:

	list of the form [dict, dict]

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: c = G.add(G.max_stable(),G.identity())
	    sage: c
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 1}
	    sage: G.unstable(c)
	    [1, 2, 3, 4]
	    sage: G.stabilize_one_step(c)
	    [{1: 2, 2: 2, 3: 3, 4: 1, 5: 3}, {1: 1, 2: 1, 3: 1, 4: 1, 5: 0}]
        """
        if config is None:
            config = self._config
        firing_vector = {}
        self.reset_config(firing_vector)
        unstable = self.unstable(config)
        for v in unstable:
            config[v] -= self.out_degree(v)
            firing_vector[v] += 1
            for e in self.outgoing_edges(v):
                if e[1] != self._sink:
                    config[e[1]]+=e[2]
        self._config = config
        return [config, firing_vector]

    def stabilize(self, config=None):
        r"""
	Returns the stabilized configuration and its firing vector.

	INPUT:

	``config`` (optional) - dict

	OUTPUT:

	list - ``[out_config, firing_vector]``

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: c = S.add(S.max_stable(), S.identity())
	    sage: S.stabilize(c)
	    [{1: 2, 2: 2, 3: 1, 4: 1, 5: 1}, {1: 1, 2: 5, 3: 7, 4: 1, 5: 6}]

	The Sandpiles current configuration is modified::

	    sage: S.config()
	    {1: 2, 2: 2, 3: 1, 4: 1, 5: 1}
        """
        if config is None:
            config = self._config
        firing_vector = {}
        self.reset_config(firing_vector)
        unstable = self.unstable(config)
        while unstable:
            firing_vector = self.add(firing_vector,self.stabilize_one_step(config)[1])
            unstable = self.unstable(config)
        self._config = config
        return [config, firing_vector]

    def new_stabilize(self, config=None):
        r"""
        Stabilize \code{config}, returning \code{[out_config, firing_vector]},
        where \code{out_config} is the modified configuration.
        """
        if config is None:
            config = self._config
        c, f = cython_stabilize(config, self.reduced_laplacian(),
            self.out_degree(), self.nonsink_vertices())
        self._config = c
        return [c, f]

    def add(self, c, d):
        r"""
        Returns the sum of ``c`` and ``d`` without stabilizing.  
	
        INPUT:

	``c``, ``d`` - dict (configurations or divisors)

	OUTPUT:

	dict (configuration or divisor)

	EXAMPLES:

	Adding configurations::

	    sage: S = sandlib('generic')
	    sage: m = S.max_stable()
	    sage: m
	    {1: 2, 2: 2, 3: 1, 4: 1, 5: 1}
	    sage: S.add(m, m)
	    {1: 4, 2: 4, 3: 2, 4: 2, 5: 2}
	
	Adding divisors::

	    sage: D = S.list_to_div([1,-3,2,0,4,-2])
	    sage: E = S.list_to_div([2, 3,4,1,1,-2])
	    sage: S.add(D,E)
	    {0: 3, 1: 0, 2: 6, 3: 1, 4: 5, 5: -4}
	"""
        sum = {}
        for v in c.keys():
            sum[v] = c[v] + d[v]
        self._config = deepcopy(sum)
        return sum

    def subtract(self, c, d):
        r"""
        Returns the difference,  ``c - d``.

	INPUT:

	``c``, ``d`` - dict (configurations or divisors)

	OUTPUT:

	dict (configuration or divisor)

	EXAMPLES:

	Subtracting configurations::

	    sage: S = sandlib('generic')
	    sage: m = S.max_stable()
	    sage: m
	    {1: 2, 2: 2, 3: 1, 4: 1, 5: 1}
	    sage: S.subtract(m, m)
	    {1: 0, 2: 0, 3: 0, 4: 0, 5: 0}
	
	Subtracting divisors::

	    sage: D = S.list_to_div([1,-3,2,0,4,-2])
	    sage: E = S.list_to_div([2, 3,4,1,1,-2])
	    sage: S.subtract(D,E)
            {0: -1, 1: -6, 2: -2, 3: -1, 4: 3, 5: 0}
        """
        sum = {}
        for v in c.keys():
            sum[v] = c[v] - d[v]
        self._config = deepcopy(sum)
        return sum

    def stably_add(self, config1, config2):
        r"""
	Returns the stabilization and corresponding firing vector of the sum of
	two configurations.

	INPUT:

        ``config1``, ``config2`` - dict (configurations)

	OUTPUT:

        list (``[stabilized sum, firing_vector]``)


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.stably_add(S.max_stable(), S.max_stable())
	    [{1: 2, 2: 2, 3: 0, 4: 1, 5: 1}, {1: 1, 2: 7, 3: 10, 4: 1, 5: 9}]
        """
        sum = {}
        for v in self._nonsink_vertices:
            sum[v] = config1[v] + config2[v]
        return self.stabilize(sum)

    def recurrent_difference(self, config1, config2):
        r"""
	Returns a recurrent configuration equivalent to the difference of the
	two given configurations, ``config1 - config2``.

	INPUT:

        ``config1``, ``config2`` - dict (configurations)

	OUTPUT:

        dict (configuration)

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: c = S.recurrent_difference(S.max_stable(), S.max_stable())
	    sage: c
	    {1: 2, 2: 2, 3: 1, 4: 1, 5: 0}
	    sage: c == S.identity()
	    True
        """
        return self.equivalent_recurrent(self.subtract(config1, config2))[0]

    def is_recurrent(self, config=None):
        r"""
        Returns True if ``config`` is recurrent, i.e., is an element of the
        sandpile group.

	INPUT:

	``config`` (optional) - dict (configuration)

	OUTPUT:

	boolean

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.is_recurrent(S.identity())
	    True
	    sage: c = S.list_to_config([0,0,0,0,0])
	    sage: S.is_recurrent(c)
	    False
        """
        if config is None:
            config = self._config
        if '_recurrents' in self.__dict__:
            return config in self._recurrents
        else:
            # add the burning configuration to config
            b = self.stably_add(config,self._burning_config)[0]
            for v in self._nonsink_vertices:
                if b[v] != config[v]:
                    return False
            return True

    def is_superstable(self, config=None):
        r"""
        Returns True if ``config`` is superstable, i.e., it does not permit
	multiset firings.

	INPUT:

	``config`` (optional) - dict (configuration)

	OUTPUT:

	boolean

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.is_superstable(S.identity())
	    False
	    sage: c = S.subtract(S.max_stable(), S.identity())
	    sage: S.is_superstable(c)
	    True
        """
        if config is None:
            config = deepcopy(self._config)
        return self.is_recurrent(self.dualize(config))


    def is_stable(self, config=None):
        r"""
        Returns True if ``config`` is stable.

	INPUT:

	``config`` (optional) - dict (configuration)

	OUTPUT:

	boolean

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.is_stable(S.max_stable())
	    True
	    sage: c = S.add(S.max_stable(), S.max_stable())
	    sage: S.is_stable(c)
	    False
        """
        if config is None:
            config = self._config
        for v in self._nonsink_vertices:
            if config[v] >= self.out_degree(v):
                return False
        return True

    def dualize(self, config=None):
        r"""
        Returns ``max_stable - config``.

        INPUT:

	``config`` (optional) - dict (configuration)

	OUTPUT:

	dict (configuration)

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.dualize(S.max_stable())
	    {1: 0, 2: 0, 3: 0, 4: 0, 5: 0}
        """
        if config is None:
            config = deepcopy(self._config)
        dual = {}
        for v in self._nonsink_vertices:
            dual[v] = self._max_stable[v] - config[v]
        self._config = dual
        return dual

    def _set_identity(self):
        r"""
	Computes ``_identity``, the variable holding the identity configuration
	of the sandpile group, when ``identity()`` is first called by a user.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_identity()
        """
        x = self.stably_add(self._max_stable, self._max_stable)[0]
        x = self.dualize(x)
        x = self.stably_add(self._max_stable, x)[0]
        self._identity = deepcopy(x)

    def identity(self):
        r"""
        Returns the identity configuration.

        INPUT: 
	
	None

        OUTPUT: 
	
	dict (the identity configuration)

	EXAMPLES::

	    sage: G = sandlib('generic')
	    sage: e = G.identity()
	    sage: x = G.stably_add(e, G.max_stable())
	    sage: x
	    [{1: 2, 2: 2, 3: 1, 4: 1, 5: 1}, {1: 1, 2: 5, 3: 7, 4: 1, 5: 6}]
	    sage: G.config_to_list(x[0])
	    [2, 2, 1, 1, 1]
	    sage: G.config_to_list(G.max_stable())
	    [2, 2, 1, 1, 1]
        """
        return deepcopy(self._identity)

    def _set_recurrents(self):
        """
	Computes ``_recurrents``, the variable holding the list of recurrent
	configurations, when ``recurrents()`` is first called by a user.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_recurrents()
        """
        save_config = deepcopy(self._config)
        self._recurrents = []
        active = [self._max_stable]
        while active != []:
            c = active.pop()
            self._recurrents.append(c)
            for v in self._nonsink_vertices:
                cnext = deepcopy(c)
                cnext[v] += 1
                self.stabilize(cnext)
                if (cnext not in active) and (cnext not in self._recurrents):
                    active.append(cnext)
        self._config = save_config

    def recurrents(self, verbose=True):
        r"""
	Returns the list of recurrent configurations as dictionaries if
	``verbose`` is ``True``, otherwise as lists of integers.

        INPUT: 
	
	``verbose`` (optional) -- boolean

        OUTPUT: 
	
	list (of recurrent configurations)


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.recurrents()
	    [{1: 2, 2: 2, 4: 1, 4: 1, 5: 1},
	     {1: 2, 2: 2, 3: 0, 4: 1, 5: 1},
	     {1: 0, 2: 2, 3: 1, 4: 1, 5: 0},
	     {1: 0, 2: 2, 3: 1, 4: 1, 5: 1},
	     {1: 1, 2: 2, 3: 1, 4: 1, 5: 1},
	     {1: 1, 2: 2, 3: 0, 4: 1, 5: 1},
	     {1: 2, 2: 2, 3: 1, 4: 0, 5: 1},
	     {1: 2, 2: 2, 3: 0, 4: 0, 5: 1},
	     {1: 2, 2: 2, 3: 1, 4: 0, 5: 0},
	     {1: 1, 2: 2, 3: 1, 4: 1, 5: 0},
	     {1: 1, 2: 2, 3: 1, 4: 0, 5: 0},
	     {1: 1, 2: 2, 3: 1, 4: 0, 5: 1},
	     {1: 0, 2: 2, 3: 0, 4: 1, 5: 1},
	     {1: 2, 2: 2, 3: 1, 4: 1, 5: 0},
	     {1: 1, 2: 2, 3: 0, 4: 0, 5: 1}]
	    sage: S.recurrents(false)
	    [[2, 2, 1, 1, 1],
	     [2, 2, 0, 1, 1],
	     [0, 2, 1, 1, 0],
	     [0, 2, 1, 1, 1],
	     [1, 2, 1, 1, 1],
	     [1, 2, 0, 1, 1],
	     [2, 2, 1, 0, 1],
	     [2, 2, 0, 0, 1],
	     [2, 2, 1, 0, 0],
	     [1, 2, 1, 1, 0],
	     [1, 2, 1, 0, 0],
	     [1, 2, 1, 0, 1],
	     [0, 2, 0, 1, 1],
	     [2, 2, 1, 1, 0],
	     [1, 2, 0, 0, 1]]
        """
        if verbose:
            return deepcopy(self._recurrents)
        else:
            verts = self.nonsink_vertices()
            return [[r[v] for v in verts] for r in self._recurrents]

    def _set_superstables(self):
        r"""
	Computes ``_superstables``, the variable holding the list of superstable
	configurations, when ``superstables()`` is first called by a user.  

	INPUT:

	None

	OUTPUT:
	
	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_superstables()
        """
        save_config = deepcopy(self._config)
        self._superstables = [self.dualize(c) for c in self._recurrents]
        self._config = save_config

    def superstables(self, verbose=True):
        r"""
	Returns the list of superstable configurations as dictionaries if
	``verbose`` is ``True``, otherwise as lists of integers.  The
	superstables are also known as `G`-parking functions.

        INPUT: 
	
	``verbose`` - boolean

        OUTPUT: 
	
	list (of superstable elements)

        
	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.superstables()
	    [{1: 0, 2: 0, 3: 0, 4: 0, 5: 0},
	     {1: 0, 2: 0, 3: 1, 4: 0, 5: 0},
	     {1: 2, 2: 0, 3: 0, 4: 0, 5: 1},
	     {1: 2, 2: 0, 3: 0, 4: 0, 5: 0},
	     {1: 1, 2: 0, 3: 0, 4: 0, 5: 0},
	     {1: 1, 2: 0, 3: 1, 4: 0, 5: 0},
	     {1: 0, 2: 0, 3: 0, 4: 1, 5: 0},
	     {1: 0, 2: 0, 3: 1, 4: 1, 5: 0},
	     {1: 0, 2: 0, 3: 0, 4: 1, 5: 1},
	     {1: 1, 2: 0, 3: 0, 4: 0, 5: 1},
	     {1: 1, 2: 0, 3: 0, 4: 1, 5: 1},
	     {1: 1, 2: 0, 3: 0, 4: 1, 5: 0},
	     {1: 2, 2: 0, 3: 1, 4: 0, 5: 0},
	     {1: 0, 2: 0, 3: 0, 4: 0, 5: 1},
	     {1: 1, 2: 0, 3: 1, 4: 1, 5: 0}]
	    sage: S.superstables(false)
	    [[0, 0, 0, 0, 0],
	     [0, 0, 1, 0, 0],
	     [2, 0, 0, 0, 1],
	     [2, 0, 0, 0, 0],
	     [1, 0, 0, 0, 0],
	     [1, 0, 1, 0, 0],
	     [0, 0, 0, 1, 0],
	     [0, 0, 1, 1, 0],
	     [0, 0, 0, 1, 1],
	     [1, 0, 0, 0, 1],
	     [1, 0, 0, 1, 1],
	     [1, 0, 0, 1, 0],
	     [2, 0, 1, 0, 0],
	     [0, 0, 0, 0, 1],
	     [1, 0, 1, 1, 0]]
        """
        if verbose:
            return deepcopy(self._superstables)
        else:
            verts = self.nonsink_vertices()
            return [[s[v] for v in verts] for s in self._superstables]

    def _set_elementary_divisors(self):
        r"""
        Computes the variable holding the elementary divisors of the sandpile
	group when ``elementary_divisors()`` is first called by the user.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_elementary_divisors()
        """
        # Sage seems to have issues with computing the Smith normal form and
        # elementary divisors of a sparse matrix, so we have to convert:
        A = [list(i) for i in self.reduced_laplacian()]
        B = matrix(ZZ,self.num_verts()-1,A)
        self._elementary_divisors = B.elementary_divisors()

    def elementary_divisors(self):
        r"""
        The elementary divisors of the sandpile group (a finite
        abelian group).

	INPUT:

	None

	OUTPUT:

	list of integers

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.elementary_divisors()
	    [1, 1, 1, 1, 15]
        """
        return deepcopy(self._elementary_divisors)

    def _set_hilbert_function(self):
        """
        Computes the variables holding the Hilbert function of the homogeneous
        homogeneous sandpile ideal, the first differences of the Hilbert
        function, and the postulation number for the zero-set of the sandpile
        ideal when any one of these is called by the user.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_hilbert_function()
        """
        v = [sum(i.values()) for i in self._superstables]
        self._postulation = max(v)
        self._first_diffs_hilb = [v.count(i) for i in range(self._postulation+1)]
        self._hilbert_function = [1]
        for i in range(self._postulation):
            self._hilbert_function.append(self._hilbert_function[i]
                +self._first_diffs_hilb[i+1])

    def first_diffs_hilb(self):
        r"""
        Returns the first differences of the Hilbert function of the homogeneous
        sandpile ideal.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	list of nonnegative integers


	EXAMPLES::

	    sage: S = sandlib('generic')
            sage: S.hilbert_function()
            [1, 5, 11, 15]
	    sage: S.first_diffs_hilb()
	    [1, 4, 6, 4]
        """
        return deepcopy(self._first_diffs_hilb)

    def hilbert_function(self):
        r"""
        Returns the Hilbert function of the homogeneous sandpile ideal.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	list of nonnegative integers

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.hilbert_function()
	    [1, 5, 11, 15]
        """
        return deepcopy(self._hilbert_function)

    def postulation(self):
        r"""
	Returns the postulation number of the sandpile ideal.  This is the
	largest weight of a superstable configuration of the graph.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	nonnegative integer


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.postulation()
	    3
        """
        return self._postulation

    def add_random(self, config=None):
        r"""
        Add one grain of sand to a random nonsink vertex.  The ``config``
        argument can be a configuration or a divisor.

        INPUT:

        ``config`` (optional) - dict (a configuration or divisor)

        OUTPUT:

        None

        EXAMPLES:

	We compute the 'sizes' of the avalanches caused by adding random grains
	of sand to the maximal stable configuration on a grid graph.  The
	function ``stabilize()`` returns the firing vector of the
	stabilization, a dictionary whose values say how many times each vertex
	fires in the stabilization.

        ::

	    sage: S = Sandpile(grid(4,4),'sink')
	    sage: a = [sum(S.stabilize(S.add_random())[1].values())\ 
	    for i in range(1000)]
	    sage: b = [[log(i+1),log(a.count(i))]\  
	    for i in [0..max(a)] if a.count(i)]
	    sage: list_plot(b)

        Also works for divisors::

            sage: S = sandlib('generic')
            sage: D = S.list_to_div([0,0,0,0,0,0])
            sage: S.add_random(D)  
            sage: D #random
            {0: 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 1}
        """
        if config is None:
            config = self._config
        C = CombinatorialClass()
        C.list = lambda: config.keys()
        config[C.random_element()] += 1

    def list_to_config(self, L):
        r"""
        Convert a list into a configuration.

        INPUT: 
	
	``L`` - list of nonnegative integers
	
        OUTPUT: 
	
	dict (configuration)


	EXAMPLES::
           
            sage: S = sandlib('generic')
            sage: S.list_to_config([3,6,4,2,7])
            {1: 3, 2: 6, 3: 4, 4: 2, 5: 7}

	NOTES:

        The method checks the length of ``L`` but not the types of its
        entries.	
        """
        if len(L) == self.num_verts()-1:
            config = {}
            for i in range(self.num_verts()-1):
                config[self._nonsink_vertices[i]] = L[i]
            return config
        else:
            raise SyntaxError, L

    def config_to_list(self, config):
        r"""
        Convert a configuration into a list in the order determined by
	``nonsink_vertices()``.

        INPUT: 
	
	``config`` - dict (configuration)
	
        OUTPUT: 
	
	list of integers


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: c = S.list_to_config([5,4,3,2,1])
	    sage: c
	    {1: 5, 2: 4, 3: 3, 4: 2, 5: 1}
	    sage: S.config_to_list(c)
	    [5, 4, 3, 2, 1]

	NOTES:

        The method checks the length of ``config`` but not the types of its
        entries.	
        """
        if len(config) == self.num_verts()-1:
            verts = self.nonsink_vertices()
            return [config[v] for v in verts]
        else:
            raise SyntaxError, config

    def equivalent_recurrent(self, config):
        """
	Returns the recurrent configuration equivalent to the given
	configuration and also returns the corresponding firing vector.

        INPUT: 
	
	``config`` - dict (configuration)

        OUTPUT: 
	
	``[configuration, firing_vector]``


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: c = [0,0,0,0,0]
	    sage: config = S.list_to_config([0,0,0,0,0])
	    sage: x = S.equivalent_recurrent(config)
	    sage: x[0] == S.identity()
	    True
	    sage: r = [x[0][v] for v in S.nonsink_vertices()]
	    sage: f = [x[1][v] for v in S.nonsink_vertices()]
	    sage: vector(r) == vector(c) - vector(f)*S.reduced_laplacian()
	    True

	NOTES:

	Let `L` be the reduced laplacian, `c` the initial configuration, `r` the
	returned configuration, and `f` the firing vector.  Then `r = c - f\cdot
	L`.
        """
        old = config
        firing_vector = {}
        self.reset_config(firing_vector)
        done = False
        while not done:
            firing_vector = self.subtract(firing_vector, self._burning_script)
            new, new_fire = self.stably_add(old, self._burning_config)
            firing_vector = self.add(firing_vector, new_fire)
            if new == old:
                done = True
            else:
                old = new
        return [new, firing_vector]

    def equivalent_superstable(self, config):
        """
	Returns the superstable configuration equivalent to the given
	configuration and also returns the corresponding firing vector.

        INPUT: 
	
	``config`` - dict (configuration)

        OUTPUT: 
	
	``[configuration, firing_vector]``


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: m = S.max_stable()
	    sage: x = S.equivalent_superstable(m)
	    sage: S.is_superstable(x[0])
	    True
	    sage: s = [x[0][v] for v in S.nonsink_vertices()]
	    sage: f = [x[1][v] for v in S.nonsink_vertices()]
	    sage: m = [m[v] for v in S.nonsink_vertices()]
	    sage: vector(s) == vector(m) - vector(f)*S.reduced_laplacian()


	NOTES:

	Let `L` be the reduced laplacian, `c` the initial configuration, `s` the
	returned configuration, and `f` the firing vector.  Then `s = c - f\cdot
	L`.
        """

        r, fv = self.equivalent_recurrent(self.dualize(config))
        zero = {}
        self.reset_config(zero)
        return [self.dualize(r), self.subtract(zero, fv)]

    def all_k_config(self, k):
        """
        The configuration having ``k`` grains of sand on each vertex.

        INPUT: 
	
	``k`` - integer
	
        OUTPUT: 
	
	dict (configuration)


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.all_k_config(7)
	    {1: 7, 2: 7, 3: 7, 4: 7, 5: 7}
        """
        return self.list_to_config([k]*len(self._nonsink_vertices))

    def compare_configs(self, c, d):
        r"""
        Returns True if each ``c`` is at least as large as ``d`` at each
        vertex.

        INPUT: 
	
	``c``, ``d`` - dict (configurations or divisors)

        OUTPUT: 
	
	boolean

        EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: c = S.all_k_config(5)
	    sage: d = S.all_k_config(4)
	    sage: S.compare_configs(c,d)
	    True
	    sage: S.compare_configs(d,c)
	    False
        """
        for v in c.keys():
            if c[v] < d[v]:
                return False
        return True

#################################################
############ Functions for divisors #############
#################################################


    def list_to_div(self, L):
        r"""
        Convert a list of integers into a divisor.

        INPUT: 
	
	``L`` - list of integers
	
        OUTPUT: 
	
	dict (configuration)


	EXAMPLES::
            
	    sage: S = sandlib('generic')
	    sage: S.list_to_div([6,5,4,3,2,1])
            {0: 6, 1: 5, 2: 4, 3: 3, 4: 2, 5: 1}

	NOTES:

        The method checks the length of ``L`` but not the types of its
        entries.	
        """
        if len(L) == self.num_verts():
            div = {}
            verts = self.vertices()
            for i in range(self.num_verts()):
                div[verts[i]] = L[i]
            return div
        else:
            raise SyntaxError, L


    def div_to_list(self, div):
        r"""
        Convert a divisor into a list in the order determined by
	``vertices()``.

        INPUT: 
	
	``div`` - dict (divisor)
	
        OUTPUT: 
	
	list of integers


	EXAMPLES::
          
	    sage: S = sandlib('generic')
	    sage: c = S.list_to_div([6,5,4,3,2,1])
	    sage: c
            {0: 6, 1: 5, 2: 4, 3: 3, 4: 2, 5: 1}
	    sage: S.div_to_list(c)
	    [6, 5, 4, 3, 2, 1]

	NOTES:

        The method checks the length of ``div`` but not the types of its
        entries.	
        """
        if len(div) == self.num_verts():
            verts = self.vertices()
            return [div[v] for v in verts]
        else:
            raise SyntaxError, div

    def all_k_div(self, k):
        """
        The divisor having ``k`` grains of sand on each vertex.

        INPUT: 
	
	``k`` - integer
	
        OUTPUT: 
	
	dict (divisor)


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.all_k_div(7)
            {0: 7, 1: 7, 2: 7, 3: 7, 4: 7, 5: 7}
        """
        return self.list_to_div([k]*len(self.vertices()))

    def linear_system(self, div):
        r"""
	Returns the complete linear system of a divisor.
	
        INPUT:

        ``div`` - divisor

        OUTPUT:
	
        dict - ``{num_homog: int, homog:list, num_inhomog:int, inhomog:list}``

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: D = S.list_to_div([0,0,0,0,0,2])
	    sage: S.linear_system(D)              
	    {'homog': [[1, 0, 0, 0, 0, 0], [-1, 0, 0, 0, 0, 0]],
	     'inhomog': [[0, 0, -1, -1, 0, -2], [0, 0, 0, 0, 0, -1], [0, 0, 0, 0, 0, 0]],
	     'num_homog': 2,
	     'num_inhomog': 3}

        NOTES:

	If `L` is the Laplacian, an arbitrary `v` such that `v\cdot L\geq -D`
	has the form `v = w + t` where `w` is in ``inhomg`` and `t` is in the
	integer span of ``homog`` in the output of ``linear_system(D)``.  

	WARNING:

	This method requires 4ti2.  After local installation of 4ti2, set the
	``path_to_zsolve`` at the beginning of ``sandpile.sage``.
        """
        import os      # should this be here?  Is there a problem with importing
        import string  # twice?  Is this even necessary with sage?

        L = self._laplacian.transpose()
        n = self.num_verts()

        # temporary file names
        lin_sys = tmp_filename()
        lin_sys_mat = lin_sys + '.mat'
        lin_sys_rel = lin_sys + '.rel'
        lin_sys_rhs = lin_sys + '.rhs'
        lin_sys_sign= lin_sys + '.sign'
        lin_sys_zhom= lin_sys + '.zhom'
        lin_sys_zinhom= lin_sys + '.zinhom'
        lin_sys_log = lin_sys + '.log'

        mat_file = open(lin_sys_mat,'w')
        mat_file.write(str(n)+' ')
        mat_file.write(str(n)+'\n')
        for r in L:
            mat_file.write(string.join(map(str,r)))
            mat_file.write('\n')
        mat_file.close()
        # relations file
        rel_file = open(lin_sys_rel,'w')
        rel_file.write('1 ')
        rel_file.write(str(n)+'\n')
        rel_file.write(string.join(['>']*n))
        rel_file.write('\n')
        rel_file.close()
        # right-hand side file
        rhs_file = open(lin_sys_rhs,'w')
        rhs_file.write('1 ')
        rhs_file.write(str(n)+'\n')
        rhs_file.write(string.join([str(-i) for i in self.div_to_list(div)]))
        rhs_file.write('\n')
        rhs_file.close()
        # sign file
        sign_file = open(lin_sys_sign,'w')
        sign_file.write('1 ')
        sign_file.write(str(n)+'\n')
        """
	Conjecture: taking only 1s just below is OK, i.e., looking for solutions
	with nonnegative entries.  The Laplacian has kernel of dimension 1,
	generated by a nonnegative vector.  I would like to say that translating
	by this vector, we transform any solution into a nonnegative solution.
	What if the vector in the kernel does not have full support though?
        """
        sign_file.write(string.join(['2']*n))  # so maybe a 1 could go here
        sign_file.write('\n')
        sign_file.close()
        # compute
        try:
            os.system(path_to_zsolve+'zsolve -q ' + lin_sys + ' > ' + lin_sys_log)
            # process the results
            zhom_file = open(lin_sys_zhom,'r')
        except IOError:
            print """
                 **********************************
                 *** This method requires 4ti2. ***
                 **********************************

            To make 4ti2 usable from Sage Sandpiles:

            * download the program from the 4ti2 homepage (http://www.4ti2.de/), and
              follow the installation instructions given there
            * open sandpiles.sage in your favorite text editor and edit the following
              line (near the beginning of the file, near the copyright statement and the
              start of the definition of the Sandpile class), replacing the
              path_to_zsolve
              string with the path to the executables in your 5ti2 directory:

                  # set the following if 4ti2 is installed
                  path_to_zsolve = '/home/davidp/math/sandpile/4ti2/linux_x86/'

            * start Sage and load sandpiles.sage.
            """
            return
        ## first, the cone generators (the homogeneous points)
        a = zhom_file.read()
        zhom_file.close()
        a = a.split('\n')
        # a starts with two numbers. We are interested in the first one
        num_homog = int(a[0].split()[0])
        homog = [map(int,i.split()) for i in a[2:-1]]
        ## second, the inhomogeneous points
        zinhom_file = open(lin_sys_zinhom,'r')
        b = zinhom_file.read()
        zinhom_file.close()
        b = b.split('\n')
        num_inhomog = int(b[0].split()[0])
        inhomog = [map(int,i.split()) for i in b[2:-1]]
        return {'num_homog':num_homog, 'homog':homog,
                'num_inhomog':num_inhomog, 'inhomog':inhomog}

    def effective_div(self, D):
        r"""
        Computes all of the effective divisors linearly equivalent to ``D``.

        INPUT:

        ``D`` - dict (divisor)

        OUTPUT:

        list (of divisors)

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: D = S.list_to_div([0,0,0,0,0,2])
	    sage: S.effective_div(D)              
	    [{0: 1, 1: 0, 2: 0, 3: 1, 4: 0, 5: 0},
	     {0: 0, 1: 0, 2: 1, 3: 1, 4: 0, 5: 0},
	     {0: 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 2}]
	    sage: [S.div_to_list(d) for d in _]
	    [[1, 0, 0, 1, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 0, 2]]
        """
        result = []
        r = self.linear_system(D)
        d = vector(self.div_to_list(D))
        for x in r['inhomog']:
            c = vector(x)*self._laplacian + d
            c = self.list_to_div(list(c))
            if c not in result:
                result.append(c)
        return result

    def r_of_D(self, D):
        r"""
        Computes ``r(D)`` and an effective divisor ``F`` such
        that ``|D - F|`` is empty.

        INPUT:

        ``D`` - divisor

        OUTPUT:

        tuple ((integer ``r(D)``, divisor ``F``))

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: D = S.list_to_div([0,0,0,0,0,4])
	    sage: E = S.r_of_D(D)                 
	    0
	    1
	    sage: E                               
	    (1, {0: 0, 1: 1, 2: 0, 3: 1, 4: 0, 5: 0})
	    sage: F = E[1]                        
	    sage: S.div_to_list(S.subtract(D,F))        
	    [0, -1, 0, -1, 0, 4]
	    sage: S.effective_div(S.subtract(D,F))
	    []
	    sage: S.r_of_D(S.list_to_div([0,0,0,0,0,-4]))
	    (-1, {0: 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: -4})
        """
        eff = self.effective_div(D)
        n = self.num_verts()
        r = -1
        if eff == []:
            return r, D
        else:
            d = vector(self.div_to_list(D))
            # standard basis vectors
            e = []
            for i in range(n):
                v = vector([0]*n)
                v[i] += 1
                e.append(v)
            level = [vector([0]*n)]
            while True:
                r += 1
                print r
                new_level = []
                for v in level:
                    for i in range(n):
                        w = v + e[i]
                        if w not in new_level:
                            new_level.append(w)
                            C = d - w
                            C = self.list_to_div(list(C))
                            eff = self.effective_div(C)
                            if eff == []:
                                return r, self.list_to_div(list(w))
                level = new_level

    def support(self, D):
        r"""
        The input is a dictionary of integers.  The output is a list of keys
        of nonzero values of the dictionary.

        INPUT: 
	
	``D`` - dict (configuration or divisor)

        OUTPUT: 
	
	list - support of the configuration or divisor


	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: c = S.identity()
	    sage: S.config_to_list(c)
	    [2, 2, 1, 1, 0]
	    sage: S.support(c)
	    [1, 2, 3, 4]
	    sage: S.vertices()
	    [0, 1, 2, 3, 4, 5]
        """
        return [i for i in D.keys() if D[i] !=0]

    def Dcomplex(self, D):
        r"""
        Returns the simplicial complex determined by the supports of the
        effective divisors linearly equivalent to ``D``.

        INPUT: 

        ``D`` - dict (divisor)

        OUTPUT: 

        simplicial complex

        EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: p = S.Dcomplex(S.list_to_div([0,1,2,0,0,1]))
	    sage: p.homology()
	    {0: 0, 1: Z x Z, 2: 0, 3: 0}
	    sage: p.f_vector()
	    [1, 6, 15, 9, 1]
	    sage: p.betti()
	    {0: 0, 1: 2, 2: 0, 3: 0}
        """
        simp = []
        for E in self.effective_div(D):
            supp_E = self.support(E)
            test = True
            for s in simp:
                if set(supp_E).issubset(set(s)):
                    test = False
                    break
            if test:
                simp.append(supp_E)
        result = []
        simp.reverse()
        while simp != []:
            supp = simp.pop()
            test = True
            for s in simp:
                if set(supp).issubset(set(s)):
                    test = False
                    break
            if test:
                result.append(supp)
        return SimplicialComplex(self.vertices(),result)

    def betti_complexes(self):
        r"""
	Returns a list of all the divisors with nonempty linear systems whose
	corresponding simplicial complexes have nonzero homology in some
	dimension. Each such divisors is returned with its corresponding
	simplicial complex.

        INPUT:

        None

        OUTPUT: 

        list (of pairs [divisors, corresponding simplicial complex])


        EXAMPLES::

	    sage: S = Sandpile({0:{},1:{0: 1, 2: 1, 3: 4},2:{3: 5},3:{1: 1, 2: 1}},0) 
	    sage: p = S.betti_complexes()
	    sage: p[0]
	    [{0: -8, 1: 5, 2: 4, 3: 1},
	     Simplicial complex with vertex set (0, 1, 2, 3) and facets {(1, 2), (3,)}]
	    sage: S.resolution()
	    'R <-- R^5 <-- R^5 <-- R^1'
	    sage: S.betti()
		       0     1     2     3
	    ------------------------------
		0:     1     -     -     -
		1:     -     5     5     -
		2:     -     -     -     1
	    ------------------------------
	    total:     1     5     5     1
	    sage: len(p)
	    11
        """
        results = []
        verts = self.vertices()
        r = self.recurrents()
        for D in r:
            D[self.sink()] = -sum(D.values())
            test = True
            while test:
                D[self.sink()] += 1 
                complex = self.Dcomplex(D)
                if sum(complex.betti().values()) > 0:
                    results.append([deepcopy(D), complex])
                if len(complex.maximal_faces()) == 1 and list(complex.maximal_faces()[0]) == verts:
                    test = False
        return results
 
#######################################
######### Algebraic Geometry ##########
#######################################

    def _compare_vertices(self, v, w):
        r"""
        Compare vertices based on their distance from the sink.

	INPUT:

	``v``, ``w`` - vertices

	OUTPUT:

	integer

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.vertices()
	    [0, 1, 2, 3, 4, 5]
	    sage: S.distance(3,S.sink())
	    2
	    sage: S.distance(1,S.sink())
	    1
	    sage: S._compare_vertices(1,3)
	    -1
        """
        return self.distance(v, self._sink) - self.distance(w, self._sink)

    def _set_ring(self):
        r"""
        Set up polynomial ring for the sandpile. 

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_ring()
        """
        # first order the vertices according to their distance from the sink
        verts = self.vertices()
        verts = sorted(verts, self._compare_vertices)
        verts.reverse()
    
        # variable i refers to the i-th vertex in self.vertices()
        names = [self.vertices().index(v) for v in verts]
        
        vars = '('
	for i in names:
            vars += 'x_' + str(i) + ','
        vars = vars[:-1] + ')' 
        # create the ring
        self._ring = singular.ring(0, vars, 'dp')
        self._variables = vars
        # create the ideal
        gens = []
        for i in self.nonsink_vertices():
            new_gen = 'x_' + str(self.vertices().index(i))
            new_gen += '^' + str(self.out_degree(i))
            new_gen += '-'
            for j in self._dict[i]:
                new_gen += 'x_' + str(self.vertices().index(j))
		new_gen += '^' + str(self._dict[i][j]) + '*'
            new_gen = new_gen[:-1]
            gens.append(new_gen)
        self._unsaturated_ideal = singular.ideal(gens)
       
    def _set_ideal(self):
        r"""
        Create the saturated lattice ideal for the sandpile.

        INPUT:
        
	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_ideal()
        """
        singular.setring(self._ring)
        self._ideal = self._unsaturated_ideal.sat(singular.ideal(
                     self._variables.replace(',','*')[1:-1]))[1]

    def unsaturated_ideal(self):
        r"""
        The unsaturated, homogeneous sandpile ideal.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	ideal

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.unsaturated_ideal()
	    x_1^3-x_4*x_3*x_0,
	    x_2^3-x_5*x_3*x_0,
	    x_3^2-x_5*x_2,
	    x_4^2-x_3*x_1,
	    x_5^2-x_3*x_2
	    sage: S.ideal()            
	    x_2-x_0,
	    x_3^2-x_5*x_0,
	    x_5*x_3-x_0^2,
	    x_4^2-x_3*x_1,
	    x_5^2-x_3*x_0,
	    x_1^3-x_4*x_3*x_0,
	    x_4*x_1^2-x_5*x_0^2
        """
        singular.setring(self._ring)
        return self._unsaturated_ideal

    def ideal(self):
        r"""
        The saturated, homogeneous sandpile ideal.

        INPUT:
	
	None
        
	OUTPUT: 
	
	ideal

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.ideal()
	    x_2-x_0,
	    x_3^2-x_5*x_0,
	    x_5*x_3-x_0^2,
	    x_4^2-x_3*x_1,
	    x_5^2-x_3*x_0,
	    x_1^3-x_4*x_3*x_0,
	    x_4*x_1^2-x_5*x_0^2
        """
        singular.setring(self._ring)
        return self._ideal
 
    def ring(self):
        r"""
	The ring containing the homogeneous sandpile ideal.  

        INPUT: 
	
	None
        
	OUTPUT: 
	
	ring

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.ring()
	    //   characteristic : 0
	    //   number of vars : 6
	    //        block   1 : ordering dp
	    //                  : names    x_5 x_4 x_3 x_2 x_1 x_0 
	    //        block   2 : ordering C

	NOTES:

	The indeterminate `x_i` corresponds to the `i`-th vertex as listed my
	the method ``vertices``. The term-ordering is degrevlex with
	indeterminates ordered according to their distance from the sink (larger
	indeterminates are further from the sink).  
        """
        return self._ring

    def _set_resolution(self):
        r"""
        Compute the free resolution of the homogeneous sandpile ideal.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_resolution()
	"""
        singular.setring(self._ring)
        self._resolution = singular.mres(self._ideal,0)

    def resolution(self, verbose=False):
        r"""
	This function computes a minimal free resolution of the homogeneous
	sandpile ideal.  If ``verbose`` is ``True``, then all of the mappings
	are returned.  Otherwise, the resolution is summarized.

        INPUT: 
	
	``verbose`` (optional) - boolean

        OUTPUT: 
	
	free resolution of the sandpile ideal

	EXAMPLES::

	    sage: S = Sandpile({0: {}, 1: {2: 2}, 2: {0: 4, 1: 1}}, 0)
	    sage: S.resolution()
	    'R <-- R^2 <-- R^1'
	    sage: S.resolution(verbose=True)
	    [1]:
	       _[1]=x_1^2-x_2^2
	       _[2]=x_1*x_2^3-x_0^4
	    [2]:
	       _[1]=x_1*x_2^3*gen(1)-x_0^4*gen(1)-x_1^2*gen(2)+x_2^2*gen(2)
	    [3]:
	       _[1]=0
        """
        singular.setring(self._ring)
        if verbose:
            return self._resolution
        else:
            betti = [self._resolution[i].size()
              for i in range(1, self._resolution.size() - 2)]
            r = ['R^'+str(i) for i in betti]
            return 'R <-- ' + join(r,' <-- ')

    def _set_groebner(self):
        r"""
        Computes a Groebner basis for the homogeneous sandpile ideal with
        respect to the standard sandpile ordering (see ``ring``).

        INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_groebner()
	"""
        singular.setring(self._ring)
        self._groebner = singular.groebner(self._ideal)

    def groebner(self):
        r"""
        Returns a Groebner basis for the homogeneous sandpile ideal with
        respect to the standard sandpile ordering (see ``ring``).

        INPUT: 
	
	None
        
	OUTPUT: 
	
	Groebner basis

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S.groebner()
	    x_2-x_0,
	    x_3^2-x_5*x_0,
	    x_5*x_3-x_0^2,
	    x_4^2-x_3*x_1,
	    x_5^2-x_3*x_0,
	    x_1^3-x_4*x_3*x_0,
	    x_4*x_1^2-x_5*x_0^2
        """
        singular.setring(self._ring)
        return self._groebner

    def _set_betti(self):
        r"""
        Computes the Betti numbers for the free resolution of the homogeneous
	sandpile ideal.

        INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::

	    sage: S = sandlib('generic')
	    sage: S._set_betti()
	"""
        singular.setring(self._ring)
        self._betti = singular.betti(self._resolution)

    def betti(self, verbose=True):
        r"""
        Computes the Betti table for the homogeneous sandpile ideal.  If
        ``verbose`` is ``True``, it prints the standard Betti table, otherwise, 
        it returns a less formated table.

        INPUT: 
	
	``verbose`` (optional) - boolean

        OUTPUT: 
	
	Betti numbers for the sandpile


	EXAMPLE::

	    sage: S = sandlib('generic')
	    sage: S.betti()
		       0     1     2     3     4     5
	    ------------------------------------------
		0:     1     1     -     -     -     -
		1:     -     4     6     2     -     -
		2:     -     2     7     7     2     -
		3:     -     -     6    16    14     4
	    ------------------------------------------
	    total:     1     7    19    25    16     4
        """
        singular.setring(self._ring)
        if verbose:
            print singular.eval('print(betti(%s),"betti")'%self._resolution.name())
        else:
            return self._betti

    def solve(self):
        r"""
        Computes approximations of the complex affine zeros of the sandpile
        ideal. 

        INPUT: 
	
	None
        
	OUTPUT: 
	
	list of complex numbers

	EXAMPLES::

	    sage: S = Sandpile({0: {}, 1: {2: 2}, 2: {0: 4, 1: 1}}, 0)
	    sage: S.solve()
	    [[0.707107*I - 0.707107, 0.707107 - 0.707107*I],
	     [-0.707107*I - 0.707107, 0.707107*I + 0.707107],
	     [-1*I, -1*I],
	     [I, I],
	     [0.707107*I + 0.707107, -0.707107*I - 0.707107],
	     [0.707107 - 0.707107*I, 0.707107*I - 0.707107],
	     [1, 1],
	     [-1, -1]]
	    sage: len(_)
	    8
	    sage: S.group_order()
	    8

	NOTES:

        The solutions form a multiplicative group isomorphic to the sandpile
	group.  Generators for this group are given exactly by ``points()``.
        """
        singular.setring(self._ring)
        v = [singular.var(i) for i in range(1,singular.nvars(self._ring))]
        vars = '('
        for i in v:
            vars += str(i)
            vars += ','
        vars = vars[:-1]  # to get rid of the final ,
        vars += ')'
        L = singular.subst(self._ideal,
                singular.var(singular.nvars(self._ring)),1)
        R = singular.ring(0,vars,'lp')
        K = singular.fetch(self._ring,L)
        K = singular.groebner(K)
        singular.LIB('solve.lib')
        M = K.solve(5,1)
        singular.setring(M)
        sol= singular('SOL').sage_structured_str_list()
        sol = sol[0][0]
        sol = [map(eval,[j.replace('i','I') for j in k]) for k in sol]
        return sol

    def _set_points(self):
        r"""
        Compute generators for the multiplicative group of zeros of the sandpile
	ideal.

	INPUT:

	None

	OUTPUT:

	None

	EXAMPLES::
	    sage: S = sandlib('generic')
	    sage: S._set_points()
	"""
        L = self._reduced_laplacian.transpose().dense_matrix()
        n = self.num_verts()-1;
        D, U, V = L.smith_form()
        self._points = []
        one = [1]*n
        for k in range(n):
            x = [exp(2*pi*I*U[k,t]/D[k,k]) for t in range(n)]
            if x not in self._points and x != one:
                self._points.append(x)

    def points(self):
        r"""
	Returns generators for the multiplicative group of zeros of the sandpile
	ideal.

        INPUT: 
	
	None
        
	OUTPUT: 
	
	list of complex numbers

	EXAMPLES:

	The sandpile group in this example is cyclic, and hence there is a
	single generator for the group of solutions.

        ::

	    sage: S = sandlib('generic')                        
	    sage: S.points()
	    [[e^(4/5*I*pi), 1, e^(2/3*I*pi), e^(-34/15*I*pi), e^(-2/3*I*pi)]]
        """
        return self._points

#######################################
#####  Symmetrical Configurations  ####
#######################################

    def is_symmetric(self, c, orbits):
        r"""
        This function checks if `c(v)` is constant over the vertices in each
        sublist of ``orbits``.

        INPUT: 
	
	 - ``c`` - configuration
         - ``orbit`` - list of lists of vertices

        OUTPUT: 
	
	boolean

	EXAMPLES::

	    sage: S = sandlib('kite')
	    sage: S.dict()
	    {0: {},
	     1: {0: 1, 2: 1, 3: 1},
	     2: {1: 1, 3: 1, 4: 1},
	     3: {1: 1, 2: 1, 4: 1},
	     4: {2: 1, 3: 1}}
	    sage: c = S.list_to_config([1, 2, 2, 3])
	    sage: S.is_symmetric(c,[[2,3]])
	    True
        """
        for x in orbits:
            if len(set([c[v] for v in x])) > 1:
                return false
        return True


    def symmetric_recurrents(self, orbits):
        r"""
        Returns the list of symmetric recurrent configurations.

	INPUT:

	``orbits`` - list of lists partitioning the vertices

	OUTPUT:

        list of recurrent configurations

	EXAMPLES::

	    sage: S = sandlib('kite')
	    sage: S.dict()
	    {0: {},
	     1: {0: 1, 2: 1, 3: 1},
	     2: {1: 1, 3: 1, 4: 1},
	     3: {1: 1, 2: 1, 4: 1},
	     4: {2: 1, 3: 1}}
	    sage: S.symmetric_recurrents([[1],[2,3],[4]])
	    [{1: 2, 2: 2, 3: 2, 4: 1}, {1: 2, 2: 2, 3: 2, 4: 0}]
	    sage: S.recurrents()
	    [{1: 2, 2: 2, 3: 2, 4: 1},
	     {1: 2, 2: 2, 3: 2, 4: 0},
	     {1: 2, 2: 1, 3: 2, 4: 0},
	     {1: 2, 2: 2, 3: 0, 4: 1},
	     {1: 2, 2: 0, 3: 2, 4: 1},
	     {1: 2, 2: 2, 3: 1, 4: 0},
	     {1: 2, 2: 1, 3: 2, 4: 1},
	     {1: 2, 2: 2, 3: 1, 4: 1}]

	NOTES:

	The user is responsible for ensuring that the list of orbits comes from
	a group of symmetries of the underlying graph.
        """
        sym_recurrents = []
        active = [self._max_stable]
        while active != []:
            c = active.pop()
            sym_recurrents.append(c)
            for orb in orbits:
                cnext = deepcopy(c)
                for v in orb:
                    cnext[v] += 1
                self.stabilize(cnext)
                if (cnext not in active) and (cnext not in sym_recurrents):
                    active.append(cnext)
        return deepcopy(sym_recurrents)

#######################################
######### Some test graphs ############
#######################################

def sandlib(selector=None):
    r"""
    Returns the sandpile identified by ``selector``.  If no argument is
    given, a description of the sandpiles in the sandlib is printed.

    INPUT: 
    
    ``selector`` - identifier or None

    OUTPUT: 
    
    sandpile or description

    EXAMPLES::

	    sage: sandlib()
	      Sandpiles in the sandlib:
		 kite : generic undirected graphs with 5 vertices
		 generic : generic digraph with 6 vertices
		 ci1 : complete intersection, non-DAG but equivalent to a DAG
		 riemann-roch1 : directed graph with postulation 9 and 3 maximal weight superstables
		 riemann-roch2 : directed graph with a superstable not majorized by a maximal superstable
		 gor : Gorenstein but not a complete intersection

	    sage: S = sandlib('gor')
	    sage: S.resolution()
	    'R <-- R^5 <-- R^5 <-- R^1'
    """
    # The convention is for the sink to be zero.
    sandpiles = {
        'generic':{
                   'description':'generic digraph with 6 vertices',
                   'graph':{0:{},1:{0:1,3:1,4:1},2:{0:1,3:1,5:1},3:{2:1,5:1},4:{1:1,3:1},5:{2:1,3:1}}
                  },
        'kite':{
                'description':'generic undirected graphs with 5 vertices',
                'graph':{0:{}, 1:{0:1,2:1,3:1}, 2:{1:1,3:1,4:1}, 3:{1:1,2:1,4:1},
                         4:{2:1,3:1}}
               },
        'riemann-roch1':{
                         'description':'directed graph with postulation 9 and 3 maximal weight superstables',
                         'graph':{0: {1: 3, 3: 1},
                                  1: {0: 2, 2: 2, 3: 2},
                                  2: {0: 1, 1: 1},
                                  3: {0: 3, 1: 1, 2: 1}
                                 }
                        },
        'riemann-roch2':{
                          'description':'directed graph with a superstable not majorized by a maximal superstable',
                          'graph':{
                                   0: {},
                                   1: {0: 1, 2: 1},
                                   2: {0: 1, 3: 1},
                                   3: {0: 1, 1: 1, 2: 1}
                                  }
                        },
	'gor':{
               'description':'Gorenstein but not a complete intersection',
               'graph':{
                        0: {},
			1: {0:1, 2: 1, 3: 4},
			2: {3: 5},
			3: {1: 1, 2: 1}
                       }
              },
	'ci1':{
               'description':'complete intersection, non-DAG but\
 equivalent to a DAG',
                   'graph':{0:{}, 1: {2: 2}, 2: {0: 4, 1: 1}}
              }

    }
    if selector==None:
        print
        print '  Sandpiles in the sandlib:'
        for i in sandpiles:
            print '    ', i, ':', sandpiles[i]['description']
        print
    elif selector not in sandpiles.keys():
        print selector, 'is not in the sandlib.'
    else:
        return Sandpile(sandpiles[selector]['graph'], 0)

#################################################
########## Some useful functions ################
#################################################

def grid(m, n):
    """
    The mxn grid graph.  Each nonsink vertex has degree 4.
    
    INPUT:
    ``m``, ``n`` - positive integers

    OUTPUT:
    dictionary for a sandpile with sink named ``sink``.

    EXAMPLE::

	sage: grid(3,4)
	{'sink': {},
	(1, 1): {'sink': 2, (1, 2): 1, (2, 1): 1},
	(1, 2): {'sink': 1, (1, 1): 1, (1, 3): 1, (2, 2): 1},
	(1, 3): {'sink': 1, (1, 2): 1, (1, 4): 1, (2, 3): 1},
	(1, 4): {'sink': 2, (1, 3): 1, (2, 4): 1},
	(2, 1): {'sink': 1, (1, 1): 1, (2, 2): 1, (3, 1): 1},
	(2, 2): {(1, 2): 1, (2, 1): 1, (2, 3): 1, (3, 2): 1},
	(2, 3): {(1, 3): 1, (2, 2): 1, (2, 4): 1, (3, 3): 1},
	(2, 4): {'sink': 1, (1, 4): 1, (2, 3): 1, (3, 4): 1},
	(3, 1): {'sink': 2, (2, 1): 1, (3, 2): 1},
	(3, 2): {'sink': 1, (2, 2): 1, (3, 1): 1, (3, 3): 1},
	(3, 3): {'sink': 1, (2, 3): 1, (3, 2): 1, (3, 4): 1},
	(3, 4): {'sink': 2, (2, 4): 1, (3, 3): 1}}
	sage: S = Sandpile(grid(3,4),'sink')
	sage: S.group_order()
	4140081
    """
    g = {}
    # corners first
    g[(1,1)] = {(1,2):1, (2,1):1, 'sink':2}
    g[(m,1)] = {(m-1,1):1, (m,2):1, 'sink':2}
    g[(1,n)] = {(1,n-1):1, (2,n):1, 'sink':2}
    g[(m,n)] = {(m-1,n):1, (m,n-1):1, 'sink':2}
    # top edge
    for col in range(2,n):
        g[(1,col)] = {(1,col-1):1, (1,col+1):1, (2,col):1, 'sink':1}
    # left edge
    for row in range (2,m):
        g[(row,1)] = {(row-1,1):1, (row+1,1):1, (row,2):1, 'sink':1}
    # right edge
    for row in range (2,m):
        g[(row,n)] = {(row-1,n):1, (row+1,n):1, (row,n-1):1, 'sink':1}
    # bottom edge
    for col in range(2,n):
        g[(m,col)] = {(m,col-1):1, (m,col+1):1, (m-1,col):1, 'sink':1}
    # inner vertices
    for row in range(2,m):
        for col in range(2,n):
            g[(row,col)] ={(row-1,col):1, (row+1,col):1, (row,col-1):1, (row,col+1):1}
    # the sink vertex
    g['sink'] = {}
    return g


def aztec(n):
    r"""
    The aztec diamond graph.

    INPUT: 

    n - integer

    OUTPUT:

    dictionary for the aztec diamond graph

    EXAMPLES::

        sage: aztec(2)
	{'sink': {(-3/2, -1/2): 2,
		  (-3/2, 1/2): 2,
		  (-1/2, -3/2): 2,
		  (-1/2, 3/2): 2,
		  (1/2, -3/2): 2,
		  (1/2, 3/2): 2,
		  (3/2, -1/2): 2,
		  (3/2, 1/2): 2},
	 (-3/2, -1/2): {'sink': 2, (-3/2, 1/2): 1, (-1/2, -1/2): 1},
	 (-3/2, 1/2): {'sink': 2, (-3/2, -1/2): 1, (-1/2, 1/2): 1},
	 (-1/2, -3/2): {'sink': 2, (-1/2, -1/2): 1, (1/2, -3/2): 1},
	 (-1/2, -1/2): {(-3/2, -1/2): 1,
			(-1/2, -3/2): 1,
			(-1/2, 1/2): 1,
			(1/2, -1/2): 1},
	 (-1/2, 1/2): {(-3/2, 1/2): 1, (-1/2, -1/2): 1, (-1/2, 3/2): 1, (1/2, 1/2): 1},
	 (-1/2, 3/2): {'sink': 2, (-1/2, 1/2): 1, (1/2, 3/2): 1},
	 (1/2, -3/2): {'sink': 2, (-1/2, -3/2): 1, (1/2, -1/2): 1},
	 (1/2, -1/2): {(-1/2, -1/2): 1, (1/2, -3/2): 1, (1/2, 1/2): 1, (3/2, -1/2): 1},
	 (1/2, 1/2): {(-1/2, 1/2): 1, (1/2, -1/2): 1, (1/2, 3/2): 1, (3/2, 1/2): 1},
	 (1/2, 3/2): {'sink': 2, (-1/2, 3/2): 1, (1/2, 1/2): 1},
	 (3/2, -1/2): {'sink': 2, (1/2, -1/2): 1, (3/2, 1/2): 1},
	 (3/2, 1/2): {'sink': 2, (1/2, 1/2): 1, (3/2, -1/2): 1}}
	sage: Sandpile(aztec(2),'sink').group_order()
        4542720

    NOTES:

    This is the aztec diamond graph with a sink vertex added.  Boundary
    vertices have edges to the sink so that each vertex has degree 4.
    """
    aztec = {}
    for i in range(n):
        for j in range(n-i):
            aztec[(1/2+i,1/2+j)] = {}
            aztec[(-1/2-i,1/2+j)] = {}
            aztec[(1/2+i,-1/2-j)] = {}
            aztec[(-1/2-i,-1/2-j)] = {}
    non_sinks = aztec.keys()
    aztec['sink'] = {}
    for vert in non_sinks:
        weight = abs(vert[0]) + abs(vert[1])
        x = vert[0]
        y = vert[1]
        if weight < n:
            aztec[vert] = {(x+1,y):1, (x,y+1):1, (x-1,y):1, (x,y-1):1}
        else:
            if (x+1,y) in aztec.keys():
                aztec[vert][(x+1,y)] = 1
            if (x,y+1) in aztec.keys():
                aztec[vert][(x,y+1)] = 1
            if (x-1,y) in aztec.keys():
                aztec[vert][(x-1,y)] = 1
            if (x,y-1) in aztec.keys():
                aztec[vert][(x,y-1)] = 1
            if len(aztec[vert]) < 4:
                out_degree = 4 - len(aztec[vert])
                aztec[vert]['sink'] = out_degree
                aztec['sink'][vert] = out_degree
    return aztec

def random_graph(num_verts, p=1/2, directed=True, weight_max=1):
    """
    A random weighted digraph with a directed spanning tree rooted at `0`.  If
    ``directed = False``, the only difference is that if `(i,j,w)` is an edge with
    tail `i`, head `j`, and weight `w`, then `(j,i,w)` appears also.  The result
    is returned as a Sage digraph.

    INPUT:

     - ``num_verts`` - number of vertices
     - ``p`` - probability edges occur
     - ``directed`` - True if directed
     - ``weight_max`` - integer maximum for random weights

    OUTPUT:

    random graph

    EXAMPLES::

	sage: g = random_graph(6,0.2,True,3)
	sage: S = Sandpile(g,0)
	sage: S.show(edge_labels = True)
    """

    a = digraphs.RandomDirectedGN(num_verts)
    b = graphs.RandomGNP(num_verts,p)
    a.add_edges(b.edges())
    if directed:
        c = graphs.RandomGNP(num_verts,p)
        # reverse the edges of c and add them in
        a.add_edges([(j,i,None) for i,j,k in c.edges()])
    else:
        a.add_edges([(j,i,None) for i,j,k in a.edges()])
        a.add_edges([(j,i,None) for i,j,k in b.edges()])
    # now handle the weights
    for i,j,k in a.edge_iterator():
        a.set_edge_label(i,j,ZZ.random_element(weight_max)+1)
    return a

def random_DAG(num_verts, p=1/2, weight_max=1):
    r"""
    Returns a random directed acyclic graph with ``num_verts`` vertices. 
    The method starts with the sink vertex and adds vertices one at a time.
    Each vertex is connected only to only previously defined vertices, and the
    probability of each possible connection is given by the argument ``p``.
    The weight of an edge is a random integer between ``1`` and
    ``weight_max``.

    INPUT: 
    
     - ``num_verts`` - positive integer
     - ``p`` - number between `0` and `1`
     - ``weight_max`` -- integer greater than `0`

    OUTPUT: 
    
    directed acyclic graph with sink `0`

    EXAMPLES::

        sage: S = random_DAG(5, 0.3)
    """
    g = {0:{}}
    for i in [1..num_verts]:
        out_edges = {}
        while out_edges == {}:
            for j in range(i):
                if p > random():
                    out_edges[j] = randint(1,weight_max)
        g[i] = out_edges
    return g

def glue_graphs(g, h, glue_g, glue_h):
    r"""
    Glue two graphs together.

    INPUT: 
    
     - ``g``, ``h`` - dictionaries for directed multigraphs
     - ``glue_h``, ``glue_g`` - dictionaries for a vertex

    OUTPUT: 
    
    dictionary for a directed multigraph


    EXAMPLES::

	sage: x = {0: {}, 1: {0: 1}, 2: {0: 1, 1: 1}, 3: {0: 1, 1: 1, 2: 1}}
	sage: y = {0: {}, 1: {0: 2}, 2: {1: 2}, 3: {0: 1, 2: 1}}
	sage: glue_x = {1: 1, 3: 2}
	sage: glue_y = {0: 1, 1: 2, 3: 1}
	sage: z = glue_graphs(x,y,glue_x,glue_y)
	sage: z
	{0: {},
	 'x0': {0: 1, 'x1': 1, 'x3': 2, 'y1': 2, 'y3': 1},
	 'x1': {'x0': 1},
	 'x2': {'x0': 1, 'x1': 1},
	 'x3': {'x0': 1, 'x1': 1, 'x2': 1},
	 'y1': {0: 2},
	 'y2': {'y1': 2},
	 'y3': {0: 1, 'y2': 1}}
	sage: S = Sandpile(z,0)                                             
	sage: S.first_diffs_hilb()
	[1, 6, 17, 31, 41, 41, 31, 17, 6, 1]
	sage: S.resolution()
	'R <-- R^7 <-- R^21 <-- R^35 <-- R^35 <-- R^21 <-- R^7 <-- R^1'

    NOTES:

    This method makes a dictionary for a graph by combining those for
    ``g`` and ``h``.  The sink of ``g`` is replaced by a vertex that
    is connected to the vertices of ``g`` as specified by ``glue_g``
    the vertices of ``h`` as specified in ``glue_h``.  The sink of the glued
    graph is `0`.

    Both ``glue_g`` and ``glue_h`` are dictionaries with entries of the form
    ``v:w`` where ``v`` is the vertex to be connected to and ``w`` is the weight
    of the connecting edge.
    """
    # first find the sinks of g and h
    for i in g:
        if g[i] == {}:
            g_sink = i
            break
    for i in h:
        if h[i] == {}:
            h_sink = i
            break
    k = {0: {}}  # the new graph dictionary, starting with the sink
    for i in g:
        if i != g_sink:
            new_edges = {}
            for j in g[i]:
                new_edges['x'+str(j)] = g[i][j]
            k['x'+str(i)] = new_edges
    for i in h:
        if i != h_sink:
            new_edges = {}
            for j in h[i]:
                if j == h_sink:
                    new_edges[0] = h[i][j]
                else:
                    new_edges['y'+str(j)] = h[i][j]
            k['y'+str(i)] = new_edges
    # now handle the glue vertex (old g sink)
    new_edges = {}
    for i in glue_g:
        new_edges['x'+str(i)] = glue_g[i]
    for i in glue_h:
        if i == h_sink:
            new_edges[0] = glue_h[i]
        else:
            new_edges['y'+str(i)] = glue_h[i]
    k['x'+str(g_sink)] = new_edges
    return k

######### Notes ################
"""
3. dir(Sandpile) and Sandpile.__dict__.keys() return methods and variables.  the
latter doesn't list the variables that are assigned by get_attr.
6. Does 'sat' return a minimal generating set?
7. Add a 'verbose' argument to identity() and max_stable().  Others?
10. Add comment that nothing has been Cythonized yet.
"""