first commit

lindenmayer.py

@@ -0,0 +1,196 @@
+from collections import deque
+
+class IndexedGenerator(object):
+    """Add a way to get a generator item by its index"""
+    def __init__(self, generator):
+        self.generator = generator
+        self.cache = []
+
+    def __getitem__(self, index):
+        for i in xrange( index - len(self.cache) + 1 ):
+            self.cache.append( self.generator.next() )
+        return self.cache[index]
+
+
+class LindenmayerSystem(IndexedGenerator):
+    """Base virtual class for a Lindenmayer system"""
+    def __init__(self, axiom, rules, angle, heading=0):
+        self.angle = angle
+        self.heading = heading
+        self.states = deque()
+        self.actions = {
+            'F': self.forward,
+            '+': self.right,
+            '-': self.left,
+            '[': self.save,
+            ']': self.restore,
+        }
+        super(LindenmayerSystem, self).__init__(self.lindenmayer(axiom, rules))
+    
+    def lindenmayer(self, axiom, rules):
+        rules = rules.items()
+
+        # def inside the lindenmayer function, so as to use "axiom" at instanciation
+        def apply(axiom, (symbol, replacement)):
+            return axiom.replace(symbol, replacement.lower())
+
+        while True:
+            yield axiom
+            axiom = reduce(apply, rules, axiom).upper()
+
+    def forward(self):
+        raise NotImplementedError
+
+    def right(self):
+        raise NotImplementedError
+
+    def left(self):
+        raise NotImplementedError
+
+    def save(self):
+        raise NotImplementedError
+
+    def restore(self):
+        raise NotImplementedError
+
+
+class TurtleLSystem(LindenmayerSystem):
+    """Draw a L-System using the Turtle module"""
+    def __init__(self, turtle, axiom, rules, angle, heading=0, size=1):
+        self.turtle = turtle
+        self.size = size
+        super(TurtleLSystem, self).__init__( axiom, rules, angle, heading )
+
+    def draw(self, depth):        
+        self.turtle.setheading(self.heading)
+
+        for char in self[depth]:
+            if char in self.actions:
+                self.actions[char]()
+
+    def forward(self):
+        self.turtle.forward(self.size)
+
+    def left(self):
+        self.turtle.left(self.angle)
+
+    def right(self):
+        self.turtle.right(self.angle)
+
+    def save(self):
+        x = self.turtle.xcor()
+        y = self.turtle.ycor()
+        h = self.turtle.heading()
+        self.states.append( (x, y, h) )
+
+    def restore(self):
+        turtle.up()
+        x, y, h = self.states.pop()
+        turtle.setx(x)
+        turtle.sety(y)
+        turtle.setheading(h)
+        turtle.down()
+
+
+class DumpTurtleLSystem(TurtleLSystem):
+    """Keep the set of uniques L-System segments drawn by the Turtle"""
+    def __init__(self, turtle, axiom, rules, angle, heading=0, size=1, rounding=10):
+        # using a set avoid duplicate segments
+        self.segments = set()
+        # nb of significant digits for rounding
+        self.rounding=10
+        super(DumpTurtleLSystem, self).__init__( turtle, axiom, rules, angle, heading, size )
+
+    def forward(self):
+        """Store segment coordinates and do a forward movement"""
+        # without rounding, there may be the same node with different coordinates, 
+        # because of error propagation
+        x1 = round( self.turtle.xcor(), self.rounding )
+        y1 = round( self.turtle.ycor(), self.rounding )
+        start = ( x1, y1 )
+        super(DumpTurtleLSystem, self).forward()
+        x2 = round( self.turtle.xcor(), self.rounding )
+        y2 = round( self.turtle.ycor(), self.rounding )
+        end = ( x2, y2 )
+        self.segments.add( (start,end) )
+
+    def draw(self, depth):
+        """Call the draw function, then clean the data"""
+        super(DumpTurtleLSystem, self).draw(depth)
+        self.clean()
+
+    def clean(self):
+        """Remove segments that have duplicated clones in the reverse direction 
+           (the segments is a set, that guarantees that no other clone exists)"""
+        for segment in self.segments:
+            for start,end in segment:
+                # FIXME surely faster to catch the exception than to do two search?
+                if (end,start) in self.segments:
+                    self.segments.remove( (end,start) )
+
+    def __str__(self):
+        dump = ""
+        for segment in self.segments:
+            for coords in segment:
+                for x in coords:
+                    dump += str(x)+" "
+            dump += "\n"
+        return dump
+
+
+def plot_segments( segments ):
+
+    import matplotlib.pyplot as plot
+    from matplotlib.path import Path
+    import matplotlib.patches as patches
+    
+    fig = plot.figure()
+    ax = fig.add_subplot(111)
+
+    for segment in segments:
+        start,end = segment
+        verts = [start,end,(0,0)]
+        codes = [Path.MOVETO,Path.LINETO,Path.STOP]
+        path = Path(verts, codes)
+        patch = patches.PathPatch(path, facecolor='none', lw=1)
+        ax.add_patch(patch)
+
+    ax.set_xlim(-50,50)
+    ax.set_ylim(-50,50)
+
+    plot.show()
+
+
+
+if __name__=="__main__":
+    import sys
+    
+    depth = 1
+    if len(sys.argv) > 1:
+        depth = int( sys.argv[1] )
+    
+    segment_size = 10
+    float_rounding = 10
+
+    from turtle import Turtle
+    turtle = Turtle()
+    turtle.speed('fastest')
+    penrose = DumpTurtleLSystem(turtle, 
+            axiom="[X]++[X]++[X]++[X]++[X]", 
+            rules={
+                'F': "",
+                'W': "YF++ZF----XF[-YF----WF]++",
+                'X': "+YF--ZF[---WF--XF]+",
+                'Y': "-WF++XF[+++YF++ZF]-",
+                'Z': "--YF++++WF[+ZF++++XF]--XF"
+            }, 
+            angle=36, heading=0, size=segment_size, rounding=float_rounding )
+
+    penrose.draw( depth )
+    #print penrose
+
+    #plot_segments( penrose.segments )
+
+    import tsplib
+    
+    tsplib.write_segments( penrose.segments, segment_size, depth, float_rounding, fd=sys.stdout )

tsplib.py

@@ -0,0 +1,249 @@
+import sys
+import scipy
+
+def write_segments( segments, size, depth, rounding, fd = sys.stdout,
+        node_coord_section=False, 
+        edge_data_section=False, 
+        edge_weight_section=True, 
+        display_data_section=True ):
+
+    # construct a {coords:id} dictionary
+    nodes = {}
+    nb = 0
+    for segment in segments:
+        for coords in segment:
+            if not nodes.has_key(coords):
+                nodes[coords] = nb
+                nb += 1
+
+    fd.write( "NAME : penrose3_%i\n" % depth)
+    fd.write("COMMENT : Rhombus Penrose tiling (type P3) as generated by a L-system, at depth %i\n" % depth)
+    fd.write("TYPE : TSP\n")
+    fd.write("DIMENSION : %i\n" % nb )
+    if edge_weight_section:
+        fd.write("EDGE_WEIGHT_TYPE : EXPLICIT\n")
+        fd.write("EDGE_WEIGHT_FORMAT : FULL_MATRIX\n")
+    if edge_data_section:
+        fd.write("EDGE_DATA_FORMAT : ADJ_LIST\n") # via the weight matrix?
+    if node_coord_section:
+        fd.write("NODE_COORD_TYPE : TWOD_COORDS\n") # do not work with concord
+    if display_data_section:
+        fd.write("DISPLAY_DATA_TYPE : TWOD_DISPLAY\n")
+
+    if node_coord_section:
+        fd.write("NODE_COORD_SECTION\n")
+        fmt = "%"+str(len(str(nb)))+"i %"+str(rounding)+"f %"+str(rounding)+"f\n"
+        for x,y in nodes:
+            fd.write(fmt % (nodes[(x,y)],x,y))
+
+    if edge_data_section:
+        fd.write("EDGE_DATA_SECTION\n")
+        for segment in segments:
+            start,end = segment
+            fd.write( str(nodes[start])+" "+str(nodes[end])+"\n" )
+
+    if edge_weight_section:
+        fd.write("EDGE_WEIGHT_SECTION\n")
+        # fill the weights matrix with size where necessary
+        weights = scipy.zeros((nb,nb), type(size))
+        for segment in segments:
+            start,end = segment
+            weights[nodes[start],nodes[end]] = size
+            weights[nodes[end],nodes[start]] = size
+            #fd.write(nodes[start],nodes[end]
+
+        fmt = "%"+str(len(str(size)))+"i "
+        for i in xrange(weights.shape[0]):
+            # full matrix
+            for j in xrange(weights.shape[1]):
+                fd.write(fmt % weights[i,j])
+            fd.write('\n')
+
+    if display_data_section:
+        fd.write("DISPLAY_DATA_SECTION\n")
+        fmt = "%"+str(len(str(nb)))+"i %"+str(rounding)+"f %"+str(rounding)+"f\n"
+        for x,y in nodes:
+            fd.write(fmt % (nodes[(x,y)],x,y))
+
+    fd.write("EOF\n")
+
+
+def read_tour_index( fd ):
+    tour = []
+    nb = int(fd.readline().strip())
+    for line in fd:
+        tour += line.split()
+
+    return map(int, tour)
+
+
+def read_nodes( fd ):
+    """Parse a .tsp file and returns a dictionary of nodes, of the form {id:(x,y)}"""
+    nodes = {}
+    data_section = False
+    for line in fd:
+        if line.strip() == "DISPLAY_DATA_SECTION":
+            data_section = True
+            continue
+
+        data = line.strip().split()
+        if len(data) != 3 or ( len(data) > 1 and data[0].isalpha() ):
+            data_section = False
+
+        if data_section == True:
+            nodes[ int(data[0]) ] = ( float(data[1]),float(data[2]) )
+
+    return nodes
+
+
+def read_vertices( fd ):
+    vertices = set()
+
+    data_section = False
+    i=-1
+    for line in fd:
+        if line.strip() == "EDGE_WEIGHT_SECTION":
+            data_section = True
+            i=0
+            continue
+
+        data = line.strip().split()
+        if len(data)==0 or ( len(data) >= 1 and data[0][0].isalpha() ):
+            data_section = False
+
+        if data_section == True:
+            for j in xrange(len(data)):
+                if float(data[j]) != 0 and (j,i) not in vertices:
+                    vertices.add( (i,j) )
+            i += 1
+
+    return vertices
+                
+
+def plot_segments( segments ):
+
+    import matplotlib.pyplot as plot
+    from matplotlib.path import Path
+    import matplotlib.patches as patches
+    
+    fig = plot.figure()
+    ax = fig.add_subplot(111)
+
+    for segment in segments:
+        start,end = segment
+        verts = [start,end,(0,0)]
+        codes = [Path.MOVETO,Path.LINETO,Path.STOP]
+        path = Path(verts, codes)
+        patch = patches.PathPatch(path, facecolor='none', lw=1)
+        ax.add_patch(patch)
+
+    ax.set_xlim(-50,50)
+    ax.set_ylim(-50,50)
+
+    plot.show()
+
+
+def plot_segments_tour( segments_1, segments_2 ):
+
+    import matplotlib.pyplot as plot
+    from matplotlib.path import Path
+    import matplotlib.patches as patches
+    
+    fig = plot.figure()
+    ax = fig.add_subplot(111)
+
+    for segment in segments_1:
+        start,end = segment
+        verts = [start,end,(0,0)]
+        codes = [Path.MOVETO,Path.LINETO,Path.STOP]
+        path = Path(verts, codes)
+        patch = patches.PathPatch(path, facecolor='0.5', lw=1)
+        ax.add_patch(patch)
+
+    for segment in segments_2:
+        start,end = segment
+        verts = [start,end,(0,0)]
+        codes = [Path.MOVETO,Path.LINETO,Path.STOP]
+        path = Path(verts, codes)
+        patch = patches.PathPatch(path, facecolor='1.0', lw=3)
+        ax.add_patch(patch)
+
+
+    ax.set_xlim(-50,50)
+    ax.set_ylim(-50,50)
+
+    plot.show()
+
+
+
+
+if __name__=="__main__":
+    import sys
+
+#    segments = [
+#            ( (0,0),(0,2) ),
+#            ( (0,2),(2,2) ),
+#            ( (2,2),(2,0) ),
+#            ( (2,0),(0,0) )
+#    ]
+#
+#    filename = "test.tsp"
+#    with open(filename,"w") as fd:
+#        write_segments( segments, fd=fd, size=1, depth=0, rounding=10 )
+#        write_segments( segments, fd=sys.stdout, size=1, depth=0, rounding=10 )
+#
+#    with open(filename,"r") as fd:
+#        nodes = read_nodes( fd )
+#
+#    print "Nodes: id (x, y)"
+#    for idx,node in nodes.items():
+#        print idx,node
+#
+#    with open(filename,"r") as fd:
+#        vertices = read_vertices( fd )
+#
+#    print "Segments: (x1,y1) (x2,y2)"
+#    segments = []
+#    for i1,i2 in vertices:
+#        print nodes[i1],nodes[i2]
+#        segments.append( (nodes[i1],nodes[i2]) )
+#
+#    plot_segments( segments )
+
+    finstance = sys.argv[1]
+    ftour = sys.argv[2]
+
+    print "Read nodes"
+    with open(finstance,"r") as fd:
+        nodes = read_nodes( fd )
+   
+    print "Read vertices"
+    with open(finstance,"r") as fd:
+        vertices = read_vertices( fd )
+
+    print "Build segments"
+    segments = []
+    for i1,i2 in vertices:
+        #print nodes[i1],nodes[i2]
+        segments.append( (nodes[i1],nodes[i2]) )
+
+#    print "Plot segments"
+#    plot_segments( segments )
+
+
+    print "Read tour"
+    with open(ftour,"r") as fd:
+        tour = read_tour_index( fd )
+
+    print "Build tour segments"
+    tour_segments = []
+    for i in xrange(0,len(tour)-1):
+        tour_segments.append( ( nodes[tour[i]],nodes[tour[i+1]] ) )
+        print tour_segments[-1]
+
+    tour_segments.append( ( nodes[tour[i+1]], nodes[tour[0]] ) )
+    print tour_segments[-1]
+
+    print "Plot tour segments"
+    plot_segments_tour( segments, tour_segments )
+