Ticket #7770: trac_7770_hanoi_tower_graph.patch

File trac_7770_hanoi_tower_graph.patch, 11.3 KB (added by rbeezer, 13 years ago)
  • sage/graphs/graph_generators.py

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1261872752 28800
    # Node ID 88549e9e6e60f194f5290b82f3bf3ed5810e6152
    # Parent  21efb0b3fc474972b5c7f617d99173536a3d79d0
    [mq]: hanoi
    
    diff -r 21efb0b3fc47 -r 88549e9e6e60 sage/graphs/graph_generators.py
    a b  
    875875        G = networkx.grid_graph(dim)
    876876        return graph.Graph(G, name="Grid Graph for %s"%dim)
    877877
     878    def HanoiTowerGraph(self, pegs, disks, labels=True, positions=True):
     879        r"""
     880        Returns the graph whose vertices are the states of the
     881        Tower of Hanoi puzzle, with edges representing legal moves between states.
     882
     883        INPUT:
     884
     885        - ``pegs`` - the number of pegs in the puzzle, 2 or greater
     886        - ``disks`` - the number of disks in the puzzle, 1 or greater
     887        - ``labels`` - default: ``True``, if ``True`` the graph contains
     888          more meaningful labels, see explanation below.  For large instances,
     889          turn off labels for much faster creation of the graph.
     890        - ``positions`` - default: ``True``, if ``True`` the graph contains
     891          layout information.  This creates a planar layout for the case
     892          of three pegs.  For large instances, turn off layout information
     893          for much faster creation of the graph.
     894
     895        OUTPUT:
     896
     897        The Tower of Hanoi puzzle has a certain number of identical pegs
     898        and a certain number of disks, each of a different radius.
     899        Initially the disks are all on a single peg, arranged
     900        in order of their radii, with the largest on the bottom.
     901
     902        The goal of the puzzle is to move the disks to any other peg,
     903        arranged in the same order.  The one constraint is that the
     904        disks resident on any one peg must always be arranged with larger
     905        radii lower down.
     906
     907        The vertices of this graph represent all the possible states
     908        of this puzzle.  Each state of the puzzle is a tuple with length
     909        equal to the number of disks, ordered by largest disk first.
     910        The entry of the tuple is the peg where that disk resides.
     911        Since disks on a given peg must go down in size as we go
     912        up the peg, this totally describes the state of the puzzle.
     913
     914        For example ``(2,0,0)`` means the large disk is on peg 2, the
     915        medium disk is on peg 0, and the small disk is on peg 0
     916        (and we know the small disk must be above the medium disk).
     917        We encode these tuples as integers with a base equal to
     918        the number of pegs, and low-order digits to the right.
     919
     920        Two vertices are adjacent if we can change the puzzle from
     921        one state to the other by moving a single disk.  For example,
     922        ``(2,0,0)`` is adjacent to ``(2,0,1)`` since we can move
     923        the small disk off peg 0 and onto (the empty) peg 1.
     924        So the solution to a 3-disk puzzle (with at least
     925        two pegs) can be expressed by the shortest path between
     926        ``(0,0,0)`` and ``(1,1,1)``.
     927
     928        For greatest speed we create graphs with integer vertices,
     929        where we encode the tuples as integers with a base equal
     930        to the number of pegs, and low-order digits to the right.
     931        So for example, in a 3-peg puzzle with 5 disks, the
     932        state ``(1,2,0,1,1)`` is encoded as
     933        `1\ast 3^4 + 2\ast 3^3 + 0\ast 3^2 + 1\ast 3^1 + 1\ast 3^0 = 139`.
     934
     935        For smaller graphs, the labels that are the tuples are informative,
     936        but slow down creation of the graph.  Likewise computing layout
     937        information also incurs a significant speed penalty. For maximum
     938        speed, turn off labels and layout and decode the
     939        vertices explicily as needed.  The
     940        :meth:`sage.rings.integer.Integer.digits`
     941        with the ``padsto`` option is a quick way to do this, though you
     942        may want to reverse the list that is output.
     943
     944        PLOTTING:
     945
     946        The layout computed when ``positions = True`` will
     947        look especially good for the three-peg case, when the graph is known
     948        to be planar.  Except for two small cases on 4 pegs, the graph is
     949        otherwise not planar, and likely there is a better way to layout
     950        the vertices.
     951
     952        EXAMPLES:
     953
     954        A classic puzzle uses 3 pegs.  We solve the 5 disk puzzle using
     955        integer labels and report the minimum number of moves required.
     956        Note that `3^5-1` is the state where all 5 disks
     957        are on peg 2. ::
     958
     959            sage: H = graphs.HanoiTowerGraph(3, 5, labels=False, positions=False)
     960            sage: H.distance(0, 3^5-1)
     961            31
     962
     963        A slightly larger instance. ::
     964
     965            sage: H = graphs.HanoiTowerGraph(4, 6, labels=False, positions=False)
     966            sage: H.num_verts()
     967            4096
     968            sage: H.distance(0, 4^6-1)
     969            17
     970
     971        For a small graph, labels and layout information can be useful.
     972        Here we explicitly list a solution as a list of states. ::
     973
     974            sage: H = graphs.HanoiTowerGraph(3, 3, labels=True, positions=True)
     975            sage: H.shortest_path((0,0,0), (1,1,1))
     976            [(0, 0, 0), (0, 0, 1), (0, 2, 1), (0, 2, 2), (1, 2, 2), (1, 2, 0), (1, 1, 0), (1, 1, 1)]
     977
     978        Some facts about this graph with `p` pegs and `d` disks:
     979
     980        - only automorphisms are the "obvious" ones - renumber the pegs.
     981        - chromatic number is less than or equal to `p`
     982        - independence number is `p^{d-1}`
     983
     984        ::
     985
     986            sage: H = graphs.HanoiTowerGraph(3,4,labels=False,positions=False)
     987            sage: H.automorphism_group().is_isomorphic(SymmetricGroup(3))
     988            True
     989            sage: H.chromatic_number()
     990            3
     991            sage: len(H.independent_set()) == 3^(4-1)
     992            True
     993
     994        TESTS:
     995
     996        It is an error to have just one peg (or less). ::
     997
     998            sage: graphs.HanoiTowerGraph(1, 5)
     999            Traceback (most recent call last):
     1000            ...
     1001            ValueError: Pegs for Tower of Hanoi graph should be two or greater (not 1)
     1002
     1003        It is an error to have zero disks (or less). ::
     1004
     1005            sage: graphs.HanoiTowerGraph(2, 0)
     1006            Traceback (most recent call last):
     1007            ...
     1008            ValueError: Disks for Tower of Hanoi graph should be one or greater (not 0)
     1009
     1010        .. rubric:: Citations
     1011
     1012        .. [ARETT-DOREE] Danielle Arett and Su Doree, Coloring and counting on the tower of Hanoi graphs,
     1013           Mathematics Magazine, to appear.
     1014
     1015
     1016        AUTHOR:
     1017
     1018        - Rob Beezer, (2009-12-26), with assistance from Su Doree
     1019
     1020        """
     1021
     1022        # sanitize input
     1023        from sage.rings.all import Integer
     1024        pegs = Integer(pegs)
     1025        if pegs < 2:
     1026            raise ValueError("Pegs for Tower of Hanoi graph should be two or greater (not %d)" % pegs)
     1027        disks = Integer(disks)
     1028        if disks < 1:
     1029            raise ValueError("Disks for Tower of Hanoi graph should be one or greater (not %d)" % disks)
     1030
     1031        # Each state of the puzzle is a tuple with length
     1032        # equal to the number of disks, ordered by largest disk first
     1033        # The entry of the tuple is the peg where that disk resides
     1034        # Since disks on a given peg must go down in size as we go
     1035        # up the peg, this totally describes the puzzle
     1036        # We encode these tuples as integers with a base equal to
     1037        # the number of pegs, and low-order digits to the right
     1038
     1039        # complete graph on number of pegs when just a single disk
     1040        edges = [[i,j] for i in range(pegs) for j in range(i+1,pegs)]
     1041
     1042        nverts = 1
     1043        for d in range(2, disks+1):
     1044            prevedges = edges      # remember subgraph to build from
     1045            nverts = pegs*nverts   # pegs^(d-1)
     1046            edges = []
     1047
     1048            # Take an edge, change its two states in the same way by adding
     1049            # a large disk to the bottom of the same peg in each state
     1050            # This is accomplished by adding a multiple of pegs^(d-1)
     1051            for p in range(pegs):
     1052                largedisk = p*nverts
     1053                for anedge in prevedges:
     1054                    edges.append([anedge[0]+largedisk, anedge[1]+largedisk])
     1055
     1056            # Two new states may only differ in the large disk
     1057            # being the only disk on two different pegs, thus
     1058            # otherwise being a common state with one less disk
     1059            # We construct all such pairs of new states and add as edges
     1060            from sage.combinat.subset import Subsets
     1061            for state in range(nverts):
     1062                emptypegs = range(pegs)
     1063                reduced_state = state
     1064                for i in range(d-1):
     1065                    apeg = reduced_state % pegs
     1066                    if apeg in emptypegs:
     1067                        emptypegs.remove(apeg)
     1068                    reduced_state = reduced_state//pegs
     1069                for freea, freeb in Subsets(emptypegs, 2):
     1070                    edges.append([freea*nverts+state,freeb*nverts+state])
     1071
     1072        H = graph.Graph({}, loops=False, multiedge=False)
     1073        H.add_edges(edges)
     1074
     1075
     1076        # Making labels and/or computing positions can take a long time,
     1077        # relative to just constructing the edges on integer vertices.
     1078        # We try to minimize coercion overhead, but need Sage
     1079        # Integers in order to use digits() for labels.
     1080        # Getting the digits with custom code was no faster.
     1081        # Layouts are circular (symmetric on the number of pegs)
     1082        # radiating outward to the number of disks (radius)
     1083        # Algorithm uses some combination of alternate
     1084        # clockwise/counterclockwise placements, which
     1085        # works well for three pegs (planar layout)
     1086        #
     1087        from sage.functions.trig import sin, cos, csc
     1088        if labels or positions:
     1089            mapping = {}
     1090            pos = {}
     1091            a = Integer(-1)
     1092            one = Integer(1)
     1093            if positions:
     1094                radius_multiplier = 1 + csc(pi/pegs)
     1095                sine = []; cosine = []
     1096                for i in range(pegs):
     1097                    angle = 2*i*pi/float(pegs)
     1098                    sine.append(sin(angle))
     1099                    cosine.append(cos(angle))
     1100            for i in range(pegs**disks):
     1101                a += one
     1102                state = a.digits(base=pegs, padto=disks)
     1103                if labels:
     1104                    state.reverse()
     1105                    mapping[i] = tuple(state)
     1106                    state.reverse()
     1107                if positions:
     1108                    locx = 0.0; locy = 0.0
     1109                    radius = 1.0
     1110                    parity = -1.0
     1111                    for index in range(disks):
     1112                        p = state[index]
     1113                        radius *= radius_multiplier
     1114                        parity *= -1.0
     1115                        locx_temp = cosine[p]*locx - parity*sine[p]*locy + radius*cosine[p]
     1116                        locy_temp = parity*sine[p]*locx + cosine[p]*locy - radius*parity*sine[p]
     1117                        locx = locx_temp
     1118                        locy = locy_temp
     1119                    pos[i] = (locx,locy)
     1120            # set positions, then relabel (not vice versa)
     1121            if positions:
     1122                H.set_pos(pos)
     1123            if labels:
     1124                H.relabel(mapping)
     1125
     1126        return H
     1127
    8781128    def HouseGraph(self):
    8791129        """
    8801130        Returns a house graph with 5 nodes.