Solving L-System Recursion

Here are step-by-step solutions to both the exercises presented at the end of Lindenmayer Systems.

Recursive version of a Koch curve L-system

Exercise: Modify the iterative code for a Koch curve L-system into recursive form. Make sure that your recursion preserves the original start and end points of the order 0 fractal - that is, if we have a Koch curve that begins at (-500, 0) and ends at (500, 0), then any order of the Koch curve should do the same.

Solution: The first thing to do is to start small. We can simplify the problem by putting aside the whole drawing business and just generating the string using a recursive mechanism. If we can get the right theorem, then we should be able to draw the shape without a problem.

Let’s go back to Lindenmayer’s original system, where the rules of production state A —> AB and B —> A. Here is our original iterative code:

def lSysGenerate(s, iterations):
    for i in range(iterations):
        s = lSysCompute(s)
        print('s =', s)
    return s

def lSysCompute(s):
    d = {'A': 'AB', 'B': 'A'}
    return ''.join([d.get(c) or c for c in s])

axiom = 'A'
iterations = 3
print(lSysGenerate(axiom, iterations))

>>> s = AB
>>> s = ABA
>>> s = ABAAB
>>> ABAAB

We can approach the problem intuitively and conservatively, by keeping the code that we know works and rewriting the code that no longer applies. If you examine the two functions, lSysGenerate() feeds the string s to lSysCompute() while lSysCompute() ‘does the work’. You can already see the similarity to mergeSort() and sierpinski() where we put the recursive action in one function and fed its outputs to other functions that merged and sorted lists, or drew triangles. So let’s begin by converting lSysGenerate() from an iterative to a recursive form.

To do that, let’s divide the problem into its pre- and post-recursive sections. Pre-recursively, it seems like there’s really no work that needs to be done. We don’t need to seed our frames with anything on the way to the base case. We just want to set up the right number of frames, so that we recurse the correct number of times and get the correct final theorem.

This implies that all the work should happen post-recursively. More specifically, this means that the base case returns the seed for all further computations, just as we returned [[]] for powerSet(). With the base case in hand, every pass through the original code’s for loop can be recursively restated as each successive frame’s application of the rules of production to the base case/preceding frame.

If we’ve done it right, our final returned string will have accumulated all the substitutions needed for the nth order. (I’m taking advantage of the fact that, since the recursion is linear, ‘frame’ and ‘order’ mean pretty much the same thing as ‘iterations’ did in the iterative code. Therefore I’ll substitute ‘order’ for ‘iterations’ from here on out.)

As for the base case, it’s obviously the axiom itself. In the case of the original L-system, this is A, which, along with order, we pass to lSysGenerate(). So with order as our counter, we keep decrementing until we get to the base case of order == 0. Following our standard recursive template, we can assert:

def lSysGenerate(axiom, order):
    if order == 0:
        return axiom
    else:
        return lSysGenerate(axiom, order - 1)

order = 3
print(lSysGenerate('A', order)

>>> A

Since the function doesn’t have multiple recursive calls, we know that we’re dealing with a fairly simple structure here. Nevertheless, let’s throw in some print-tracing to keep track of frames, along with splitting the recursive call from the return statement with r:

frame = 0

def lSysGenerate(axiom, order):
    global frame
    if order == 0:
        frame += 1
        print('base case, frame =', frame)
        print('returning', axiom)
        return axiom
    else:
        frame += 1
        print('frame =', frame)
        r = lSysGenerate(axiom, order - 1)
        frame -= 1
        print('frame =', frame)
        print('r =', r)
        return r

print(lSysGenerate('A', 3)

>>> frame = 1
>>> frame = 2
>>> frame = 3
>>> base case, frame = 4
>>> returning A
>>> frame = 3
>>> r = A
>>> frame = 2
>>> r = A
>>> frame = 1
>>> r = A
>>> A

We next want the opportunity to apply the rules of production to the string axiom that is being passed back to us from the base case. To do this, all we need is to call lSysCompute(), and the most concise way to do it is in the return statement itself. (Calling the external function at the return statement was exactly what we did when we called the merge() function in our mergeSort() program.)

So far we have:

def lSysGenerate(axiom, order):
    if order == 0:
        return axiom
    else:
        return lSysCompute(lSysGenerate(axiom, order - 1))

def lSysCompute(s):
    d = {'A': 'AB', 'B': 'A'}
    return ''.join([d.get(c) or c for c in s])

axiom = 'A'
order = 3
print(lSysGenerate(axiom, order))

>>> ABAAB

Ok, then. What about drawing? The fact is that we cannot draw the shape until we have the final string in hand, and this doesn’t happen until the recursive has completed. Since the iterative code doesn’t draw anything until the final theorem is generated, there is no change in the reationship between the draw() and lSysGenerate() functions.

Implementing absolute distance

But we are still missing a crucial part of the solution, which is the ability to draw at a scale that preserves the absolute distance described by the 0th order. With koch() and sierpinski(), we linked the variable tracking the fractal’s order with the length of the line that would be drawn. The higher the order, the shorter the line (or the smaller the shape).

In the current case, we can’t do anything at the base case except retrieve the axiom, and we can’t draw the complete figure until we’ve exited the recursion. So it looks like we’re asking two separate things of our code: figure out the smallest line length we need, and compute the final iteration of the L-system string. We already know (and pretty much don’t have a choice) about how the latter works. The trick is to figure out where to compute the former, and how to extract it from the recursive mechanism.

If we want two values from our recursion, it makes sense to include more than just the string in the recursive function’s return statement. This means that we will also have two variables in the base case’s return statement. Since we’ll have to return two values through the recursive cases, we’ll have to be careful about what we’re subjecting to computation. Recall from our expansion of Pascal’s triangle, we went from wanting to pass a single list (ie, the nth layer of the triangle) to a list of lists (ie, all layers of the triangle up to and including the nth layer). The trick is that, while we wanted the entire list of lists, we only wanted to recursively operate on the last returned sublist. We’ll be doing something similar here.

So far, we’ve got the string-generation part down, and there’s no need to mess with that. To find the smallest line length length, we know that for every order, length will be recomputed successively as length /= 3 . So it makes sense to conduct the division pre-recursively. Once we get down to the base case, we’ll have the correct minimum value for length. Now all we need to do is pass that value of length - unchanged - back through the recursive process and hand it to the frame that initially called the recursive function.

Here’s the final code for lSysGenerate(), which I’ll walk through below:

def lSysGenerate(axiom, order, length):
    if order == 0:
        return [axiom, length]
    else:
        length /= 3
        r = lSysGenerate(axiom, order - 1, length)
        return [lSysCompute(r[0]), r[1]]

We first need to include length as an argument when defining (and calling) lSysGenerate(). Assuming that order > 0, we skip the if block and enter the else block. There, we modify the pre-recursive portion of the else block to compute length /= 3. This gets us to the point where we recursively call the function.

Now, the point of the recursive call is to get us to the base case, with the additional requirement that we bring the successively subdivided variable length along with us. Since our function currently has the parameters lSysGenerate(axiom, order, length), we decrement order - 1 and the latest value of length travels with it.

For the base case order == 0, we don’t need to recompute length, but we do need to capture it as part of our return statement. So we’ll recast our return statement as a list, [axiom, length]. In the case of the Koch curve, if we wanted order 2 and had an initial length of 1000, what should be returned is ['A', 111].

But returned as what? Here I’ve been a bit more explicit, as I want to emphasize the fact that we’re passing a list back to the calling frame. As I’ve done throughout this guide, I use the variable name r as a simple placeholder, which gives us a more intuitive view of the return statement. Since r is a list that consists of the string and the minimum length, we want to apply the rules of production (by calling lSysCompute()) to the first index item but leave the second one untouched.

In the interest of compact code, you could also write:

return [lSysCompute(lSysGenerate(axiom, order - 1, length)[0]), lSysGenerate(axiom, order - 1, length)[1]]

But this is difficult to read, and also implies an additional and unnecessary computation. Be nice to other people, and make your code easy to read.

The last modification that we need to make is in the arguments of the draw() function, since we are now accessing two index items from a list - r[0] == theorem and r[1] == length. Here is the complete code:

import string
import turtle

def lSysGenerate(axiom, order, length):
    if order == 0:
        return [axiom, length]
    else:
        length //= 3
        r = lSysGenerate(axiom, order - 1, length)
        return [lSysCompute(r[0]), r[1]]

def lSysCompute(theorem):
    d = {'A': 'A-A++A-A'}
    return ''.join([d.get(c) or c for c in lString])

def draw(t, theorem, length, angle):
    for c in theorem:
        if c in string.ascii_letters:
            t.fd(length)
        elif c == '-':
            t.left(angle)
        elif c == '+':
            t.right(angle)

def main():
    t = turtle.Turtle()
    wn = turtle.Screen()
    wn.bgcolor('black')

    t.color('orange')
    t.pensize(1)
    t.penup()
    t.setpos(-500, 0)
    t.pendown()
    t.speed(0)

    axiom = 'A'
    length = 1000
    angle = 60
    order = 4

    r = lSysGenerate(axiom, order, length)
    draw(t, r[0], r[1], angle)

    wn.exitonclick()

main()

You can quickly test this by writing a little for loop outside of main() that overlaps each order of the curve on top of the last, using a different pen color:

def main(order):
    t = turtle.Turtle()
    wn = turtle.Screen()
    wn.bgcolor('black')

    colormap = ['red', 'orange', 'yellow', 'green', 'blue', 'violet', 'white']

    t.color(colormap[i])
    t.pensize(3)
    t.penup()
    t.setpos(-1250, -600)
    t.pendown()
    t.speed(0)

    axiom = 'A'
    length = 1000
    angle = 60

    r = lSysGenerate(axiom, order, length)
    draw(t, r[0], r[1], angle)

order = 7
for i in range(order):
    main(i)

You’ll see that the code applies to any combination of alphabet, axiom and rules of production, although with varying and sometimes surprising results.

Recursively finding a target string

Exercise: For a given L-system, find if a string target exists as part of the system after a certain number of iterations. Write a recursive solution that returns the order where target appears, as well as its index position in the string as it exists in that order.

For example, say that we want to find at what point target = 'BAABAABA' first appears in Lindenmayer’s original system. After running the function the output should print target BAABAABA is at index position 6 at order 7. If it doesn’t, output should be target BAABAABA not found.

Solution: Given the generative nature of L-systems, keep in mind that target may show up in one order but, thanks to the next application of the rules of production, be substituted out at the next! By the same token, you have to write your solution so that what’s returned is the first time target is found, not the most recent.

I mention this because recursion tends to favor retrieval of values at two points: at the base case, and at the end of the entire recursive cascade. It’s more difficult to capture values during recursion. One way we did this was with mergeSort(), where we called an outside function for input that was then passed into the recursive cascade. We used this same technique for the solution to the first exercise above, where we called lSysCompute() at every return statement of lSysGenerate().

Here is a somewhat different approach, where we don’t have to invoke another function. Instead, we’ll store the values as additional parameters of the recursive function itself. Let’s start with the solution we came up with for the previous exercise. Since we’re only interested in finding a string in the theorem, we can use Lindenmayer’s original system and dispense with a drawing component.

def lSysGenerate (axiom, order):
    if order == 0:
        return axiom
    else:
        return lSysCompute(lSysGenerate(axiom, order - 1))

def lSysCompute(lString):
    d = {'A': 'AB', 'B': 'A'}
    return ''.join([d.get(c) or c for c in lString])

def main():
    axiom = 'A'
    order = 8
    print(lSysGenerate(axiom, order))

main()

>>> ABAABABAABAABABAABABAABAABABAABAABABAABABAABAABABAABABA

Great. We already know that we have to pass another argument to lSysGenerate() - the string that we’re looking for. Let’s call it target.

 def lSysGenerate (axiom, order, target):
     if order == 0:
         return axiom
     else:
         return lSysCompute(lSysGenerate(axiom, order - 1, target))

 def lSysCompute(lString):
     d = {'A': 'AB', 'B': 'A'}
     return ''.join([d.get(c) or c for c in lString])

 def main():
     axiom = 'A'
     order = 8
     target = 'BAABAABA'
     print(lSysGenerate(axiom, order, target))

We also have to unpack the else block in lSysGenerate() so that we can test for the presence of target. If it’s found, we want to record that, so we’ll need a variable hit that we’ll initially set to None. This variable is also added as a parameter to the function, so that it’s carried along with everything else.

def lSysGenerate (axiom, order, target):
    hit = None
    if order == 0:
        return [axiom, hit]
    else:
        r = lSysGenerate(axiom, order - 1, target)
        if target in r[0]:
            r[1] = r[0].find(target)
        return [lSysCompute(r[0]), r[1]]

At the base case, we return a list, initially ['A', None], so we know that r[0] is the string, and, if we find target in r[0] is True, then the index position gets recorded in the second list item, r[1]. But we’re still missing the exact order (or frame, or iteration) at which point this happens. So we’ll add another item to our base case’s returned list, order.

def lSysGenerate (axiom, order, target):
    hit = None
    if order == 0:
        return [axiom, hit, order]
    else:
        r = lSysGenerate(axiom, order - 1, target)
        if target in r[0]:
            r[1] = r[0].find(target)
            r[2] = order
        print(r)
        return [lSysCompute(r[0]), r[1], r[2]]

I snuck in a print(r) statement to see what’s going on here:

>>> ['A', None, 0]
>>> ['AB', None, 0]
>>> ['ABA', None, 0]
>>> ['ABAAB', None, 0]
>>> ['ABAABABA', None, 0]
>>> ['ABAABABAABAAB', None, 0]
>>> ['ABAABABAABAABABAABABA', 6, 7]
>>> ['ABAABABAABAABABAABABAABAABABAABAAB', 6, 8]
>>> ['ABAABABAABAABABAABABAABAABABAABAABABAABABAABAABABAABABA', 6, 8]

Hmmm, so we’re getting the correct index position for target but the value for order isn’t sticking. This is because every frame re-checks to see whether target shows up within the namespace of r[0] for that frame. We need to be able to set things up so that, once target is found, we can stop looking, and preserve both the index position and order as they were in that frame. There is no other way, since, unlike a loop, we can’t just break out of the recursion.

To do this, we set up a variable flag. Initially set to False, once target is found, we re-set flag to True and integrate it as a condition for the if block. Once flag == True then the loop is never triggered again, even if target in r[0] continues to be True for succeeding frames. Here’s our code so far:

def lSysGenerate (axiom, order, target, flag):
    hit = None
    if order == 0:
        return [axiom, hit, order, flag]
    else:
        r = lSysGenerate(axiom, order - 1, target, flag)
        if target in r[0] and r[3] == False:
            r[1] = r[0].find(target)
            r[2] = order
            r[3] = True
        print(r)
        return [lSysCompute(r[0]), r[1], r[2], r[3]]

This looks pretty good! We are, however, missing one last piece. It’s always important to consider edge cases - what if target is part of the axiom itself? As we have written the code so far, we will only return True for order == 1 even if target is present in the 0th order. So we have to add a check at the base case. Here is the final, complete code:

def lSysGenerate (axiom, order, target, flag):
    hit = None
    if order == 0:
        if target in axiom:
            hit = axiom.find(target)
            flag = True
        return [axiom, hit, order, flag]
    else:
        r = lSysGenerate(axiom, order - 1, target, flag)
        if target in r[0] and r[3] == False:
            r[1] = r[0].find(target)
            r[2] = order
            r[3] = True
        return [lSysCompute(r[0]), r[1], r[2], r[3]]

def lSysCompute(lString):
    d = {'A': 'AB', 'B': 'A'}
    return ''.join([d.get(c) or c for c in lString])

def main():
    axiom = 'A'
    order = 8
    target = 'BAABAABA'
    flag = False
    r = lSysGenerate(axiom, order, target, flag)
    if r[1] != None:
        print(``target``, target, 'is at index position', r[1], 'at order', r[2])
    else:
        print(``target``, target, 'not found')

main()

>>> target BAABAABA is at index position 6 at order 7

This code may not be as elegant or compact as many of the high-level recursive examples. On the other hand, it shows how you can retrieve values while the recursive cascade is still unfolding. Also, by gradually building up the code, I hope I made it a little less intimidating than if I’d simply introduced the final solution straight away.