L-Systems in OCaml

Question

For reference: L-system

I'm mostly looking for general style stuff here and in particular:

• if my type parameters are too crazy — I wanted two parameters so that variables and constants could have their own types.
• if I could be using some library features that I don't know about.
• if there might be some type based way to avoid the exception that I raise in do_rule. There should always be a rule for each possible variable, so the exception isn't logically necessary, but I couldn't figure out how to tell that to the compiler.

(* ************************************************************************* *
 * L-System types
 *)

(* Some elements of L-Systems can be replaced (Variables from alphabet 'a)
   and some cannot (Constants from alphabet 'b) *)
type ('a, 'b) element =
    | Var of 'a
    | Const of 'b ;;

(* Type for current state of L-System. The system's axiom is its IC state *)
type ('a, 'b) state = (('a, 'b) element) list ;;

(* Rules associate elements of the Variable type with states *)
type ('a, 'b) rule = 'a * ('a, 'b) state ;;

(* L-Systems just have to define their axiom and production rules *)
type ('a, 'b) l_system = {
    axiom : ('a, 'b) state ;
    rules : ('a, 'b) rule list ;
} ;;


(* ************************************************************************* *
 * Production
 *)

(* Run the L-System for @n generations *)
let produce (system : ('a, 'b) l_system) (n : int) : ('a, 'b) state =
    (* helper: find and apply @x's rule *)
    let rec do_rule (x : 'a) (rules : (('a, 'b) rule) list) =
        match rules with
        | [] -> 
        raise (Failure "Incomplete rule set.")
        | (y, st)::rest -> 
        if x == y then st else do_rule x rest
    in
    (* helper: find this state's successor *)
    let rec advance state rules acc : ('a, 'b) state = 
        match state with
        | [] -> 
        acc
        | (Var x)::tl ->
        advance tl rules (acc @ (do_rule x rules))
        | (Const x)::tl ->
        advance tl rules (acc @ [Const x])
    in
    (* main recursive helper *)
    let rec iterate state rules n : ('a, 'b) state =
        match n with
        | 0 -> state
        | _ -> iterate (advance state rules []) rules (n - 1)
    in iterate system.axiom system.rules n ;;


(* ************************************************************************* *
 * Examples
 *)

(* Algae *)
type two_alphabet = X | Y ;;
let algae : (two_alphabet, unit) l_system = {
    axiom = [Var X] ;
    rules = [
             (Y, [Var X]);
             (X, [Var X; Var Y])
            ]
} ;;
let bloom = produce algae 10 ;;

(* Dragon *)
type square_draw_constants = DrawForward | Left90 | Right90 ;;
let dragon : (two_alphabet, square_draw_constants) l_system = {
    axiom = [Const DrawForward; Var X] ;
    rules = [
             (X, [Var X; Const Right90; Var Y; Const DrawForward; Const Right90]);
             (Y, [Const Left90; Const DrawForward; Var X; Const Left90; Var Y])
            ]
} ;;
let a_dragon = produce dragon 10 ;;

Answer 1

This was fun to play with, so thanks for writing it and posting the question!

After experimenting with it for a while, I tried to see if I could reduce the size of produce. I don't have all the intermediate steps to show, but I made the following transformations:

• I removed the type specifiers since the compiler was able to determine the types on its own. This is a matter of taste and you could certainly put them back in. However, it made it easier for me to see the parameters which came in handy in later transformations.

• advance is only used in iterate. And do_rule is only used in advance so I moved them to be declared inside each other.

• Once the functions were nested in their enclosing function, I was able to remove some of the parameters since they could be referenced from the outer-scope.

• Your mapping of variables to their expansion is using an "association list", so I used a function in the List module to do the lookup. If the mapping isn't there, an exception is thrown so we don't need to raise an exception ourselves (although the error message won't be as nice as yours.)

• Since we're iterating through a list and accumulating a result, I replaced your advance function with an instance of List.fold_left.

After applying these transformations, the function looks like this:

let produce system =
  let rec iterate state = function
    | 0 ->
       state
    | n ->
       let subst acc = function
         | Var x ->
            acc @ (List.assoc x system.rules)
         | (Const _) as item ->
            acc @ [item]
       in
       iterate (List.fold_left subst [] state) (n - 1)
  in iterate system.axiom

Playing with the system further showed performance problems when the iterations got to 20. The reason is that produce appends each item to the end of a list, which is an O(n) operation. Using your algae alphabet as an example, 30 expansions resulted in more than 2 million elements in the list. As that result was being built, it created temporary lists from 1 up to 2 million elements! This makes your algorithm O(n^2) which is why it works fine with a few expansions, but rapidly slows down when the iterations get larger.

We know that adding an element to the beginning of the list is a constant-time operation so, if we build the list in reverse and then call List.rev to reverse the result, we change the algorithm to O(n) complexity.

I also wondered if it was necessary to build the intermediate lists. On pen and paper, it makes sense to write down the steps but maybe we could write the algorithm so that it only generates the final list. It turns out that it's easy to do this: we start with an empty list and use List.fold_left on the initial state. Anytime we find an item to expand, we call List.fold_left with the expansion and feed it our in-progress result.

This new version, produce', looks like this:

let produce' system n =
  let rec process n res = function
    | (Var x) as h ->
       if n > 0 then
         match List.assoc x system.rules with
         | nl ->
            List.fold_left (process (n - 1)) res nl
         | exception Not_found ->
            h :: res
       else
         h :: res
    | (Const _) as h ->
       h :: res in
  List.rev @@ List.fold_left (process n) [] system.axiom

I did some rough performance measurements using the following function:

let tm f a b =
  let t_start = Unix.gettimeofday () in
  let result = f a b in
  let t_end = Unix.gettimeofday () in
  (List.length result, t_end -. t_start);;

tm produce algae 1;;  (* for example *)

Running several values for n results in the following timing on my system:

 N    produce   produce'    final size

 1       2 us       2 us             1
 5       8 us       6 us            13
10     285 us      49 us           144
15      33 ms     0.5 ms         1,597
20     7.3 s      6.4 ms        17,711
25                 90 ms       196,418
30                1.2 s      2,178,309
35                 14 s     24,157,817

I didn't run the larger numbers on produce because the running times were growing quickly and I didn't want to wait.

Hope you find this helpful.