How can I improve my code for drawing a tree?

Question

This code draws a tree, but it's not in the functional style.

tree[set_, m_] := 
  Module[{θ, mz, my, temp, k, i, p1, p2, a, b, c, d, e, f},
   θ = Pi/6;
   mz = {{Cos[θ], -Sin[θ]}, {Sin[θ], Cos[θ]}}; 
   my = {{Cos[-θ], -Sin[-θ]}, {Sin[-θ], Cos[-θ]}};
   temp = set;
   p1 = {};(*store branch*)
   p2 = {};(*store leaf*)
   For[k = 1, k <= m, k++, p2 = {};
    For[i = 1, i <= Length[temp], i++, a = temp[[i, 2]];
     b = temp[[i, 1]];
     c = 1/3*a + 2/3*b;
     d = c + mz.(a - b)*(1/3);
     e = c + my.(a - b)*(1/3);
     AppendTo[p1, {Brown, Thickness[0.035/k], Line[{b, c}]}];
     AppendTo[p2, {c, a}];
     AppendTo[p2, {c, d}];
     AppendTo[p2, {c, e}];
     ];
    temp = p2;
    ];

   Show[Graphics[{p1, {RGBColor[0.1, 0.42, 0.17], Line[p2]}}, 
     Axes -> 1, PlotRange -> {{-0.2, 0.2}, {-0.2, 1.}}]]];

tree[{{{0., 0.}, {0., 1.}}}, 7]

A friend of mine tried rewrite it with Nest, but it worked slower than the original code.

(leaf = Nest[Apply[Transpose[{Table[(2 #1 + #2)/3, {3}], 
           {RotationTransform[-5 Pi/6, (2 #1 + #2)/3][#1], 
          RotationTransform[
            5 Pi/6, (2 #1 + #2)/3][#1], #2}}] &, #, {-3}] &, 
            {{{0, 0.}, {0, 1}}}, 7]; 
  Graphics[Map[{RGBColor[0.1, 0.42, 0.17], Line[#]} &, leaf, {-3}]]) // AbsoluteTiming

I think Nest is slower than loops. How do I rewrite it to make it faster?

If possible, use Fold, Nest or similar functional constructs.

Answer 1

Here's my approach. I've put the core of the algorithm into a function iterate which takes a leaf and returns a branch and 3 leaves. Each leaf or branch is defined by a start point, an end point, and your k variable which is used to set the Thickness.

The tree is constructed by Nesting the iterate function 7 times, starting with an initial leaf expression.

Finally I use Cases to extract the branches and leaves and convert them to graphics primitives.

θ = Pi/6;
mz = {{Cos[θ], -Sin[θ]}, {Sin[θ], Cos[θ]}};
my = {{Cos[-θ], -Sin[-θ]}, {Sin[-θ], Cos[-θ]}};

iterate[leaf[{b_, a_, k_}]] := Module[{c, d, e},
  c = 1/3*a + 2/3*b;
  d = c + mz.(a - b)*(1/3);
  e = c + my.(a - b)*(1/3);
  {branch[{b, c, k + 1}], leaf[{c, a, k + 1}], leaf[{c, d, k + 1}], leaf[{c, e, k + 1}]}]

iterate[b_branch] := b

tree = Nest[Flatten[iterate /@ #] &, {leaf[{{0., 0.}, {0., 1.}, 0}]}, 7];

Graphics[{
  Brown,
  Cases[tree, branch[{x_, y_, k_}] :> {Thickness[0.035/k], Line[{x, y}]}],
  RGBColor[0.1, 0.42, 0.17], 
  Cases[tree, leaf[{x_, y_, k_}] :> Line[{{x, y}}]]
  }]

Answer 2

There are now a couple of answers that show how to generate your tree in a very functional way. Really, though, your approach is not too bad so it might be worth seeing how to make minimal changes to your algorithm and speed up the result.

In terms of speed, the most serious problem with your code is the repeated use of AppendTo. The reason is that lists are implemented internally as arrays so the Append and AppendTo commands require that the entire list be moved in memory. This is very slow, particularly when building potentially long lists of arbitrary length. Here's a simple example:

list = {};
Do[AppendTo[list, k], {k, 1, 100000}] // Timing

(* Out: {29.5157, Null} *)

If you need to build a list of arbitrary length, you have several alternatives; pre-allocation and use of Sow/Reap are commonly mentioned in this context. Another undocumented but oft mentioned tool is the Bag buried in the Internal context. You can reimplement the above using a Bag like so:

bag = Internal`Bag[];
Do[Internal`StuffBag[bag, k], {k, 1, 100000}] // Timing

(* Out: {0.044028, Null} *)

It's so much faster since the internal representation is completely different - a linked list, as opposed to an array I suppose. If we're going to use the Bag commands a lot, we can add the Internal context to our $ContextPath.

If[! MemberQ[$ContextPath, "Internal`"],
    AppendTo[$ContextPath, "Internal`"]];

Now, here's a rewrite of your algorithm using Bags in place of Lists. There are a couple of other small differences as well; I've replaced the For loops with Do loops, for example, but the performance change is primarily due to the Bag/List difference.

tree2[theta_, depth_,initialDirection_:{{{0., 0.}, {0., 1.}}}] := 
  Module[{rotateLeft, rotateRight, k, i, temp, p1, p2, a, b, c, d, e},
   rotateLeft = RotationMatrix[theta];
   rotateRight = RotationMatrix[-theta];
   temp = Bag[{{{0., 0.}, {0., 1.}}}];
   p1 = Bag[];(*store branch*)
   Do[
     p2 = Bag[];(*store leaf*)
     Do[
       {b,a} = BagPart[temp,i];
       c = a/3+2b/3;
       d = c + rotateLeft.(a-b)/3;
       e = c + rotateRight.(a-b)/3;
       StuffBag[p1, {Brown, Thickness[0.035/k], Line[{b, c}]}];
       StuffBag[p2, {c,a}];
       StuffBag[p2, {c,d}];
       StuffBag[p2, {c,e}],
       {i, 1, BagLength[temp]}];
     temp = p2,
    {k, 1, depth}];
    Graphics[{BagPart[p1,All], 
     {RGBColor[0.1, 0.42, 0.17], Line[BagPart[p2,All]]}}, 
     Axes -> 1, PlotRange->{{-0.2,0.2},{-0.2,1.}}]
];

Let's try it at depth 9 while supressing the output so that we time just the computation.

treePic2 = tree2[Pi/6, 8]; // Timing

(* Out: {0.167719, Null} *)

Here's the same computation for your tree command.

treePic = tree[{{{0.0, 0.0}, {0.0, 1.0}}}, 8]; // Timing

(* Out: {7.135478, Null} *)

At depth 10, the Bag version takes about 1.5 seconds; I aborted the List version. You can see the images by just executing the name:

treePic2

(* Out: Same pic!!! *)

Answer 3

Your tree can be described using a parametric L-system, as described in Prusinkiewicz and Lindenmayer's classic The Algorithmic Beauty of Plants. Here is a very short description of a Mathematica implementation. A more complete description can be found in this paper; largely, this is an extension of the ideas in Chapter 8 of Stan Wagon's Mathematica in Action. As the paper was written with V5, it could certainly be updated, but here are the basics.

In the language of L-systems, the basic structure your tree can be described something like so:

subs = F[s_, w_] (* initial configuration, one step forward *) :> {
  F[0.6 s, 0.7 w] (* Forward, not as far and not as thick *),
  push (* store current state *),
  r[Pi/6] (* Rotate left *),
  F[0.2 s, 0.7 w] (* Another step forward *),
  pop  (* Recall last position stored *),
  (* Then more of the same *)
  push, F[0.4 s, 0.7 w], pop, push, r[-Pi/6], 
  F[0.2 s, 0.7 w], pop, 
  setWidth[w]
};

Now, we iterate this replacement rule starting from a single step forward and take a look at a small portion of the result. (Fundamentally, this is where the iteration occurs.)

instructions = Flatten[Nest[# /. subs &, {setWidth[0.035], F[1, 0.035]}, 7]];
instructions[[10 ;; 30]]

(* Out: {push, r[-Pi/6], F[0.4, 0.00288], pop, setWidth[0.00411], push, r[Pi/6],
 F[0.4, 0.00288], push, r[Pi/6], F[0.16, 0.00288], pop, push, F[0.28,0.00288], 
 pop, push, r[-Pi/6], F[0.16, 0.00288], pop, setWidth[0.00411], pop} *)

Next, we need to interpret these instructions as a list of grapics primitives.

theta = Pi/6;
thePath = {};
lastpt = {0,0};
dir = {0,1};
currentWidth = 0.04;
stack = N[{lastpt, dir, currentWidth}];
instructions /. {
  (* Movement *)
  F[x_, t___] :> 
    (thePath = {thePath, {currentWidth, {lastpt, lastpt += x*dir}}};),
  (* Rotation *)
  r[t_] :> (dir = RotationMatrix[t].dir;),
  (* Branching structures *)
  push :> (stack = Flatten[{{{lastpt, dir,
       rotleft, rotright, currentWidth, currentColor}}, 
       stack}, 1];),
  pop :> (dir = stack[[1, 2]]; lastpt = stack[[1, 1]]; 
      rotleft = stack[[1, 3]]; rotright = stack[[1, 4]]; 
      currentWidth = stack[[1,5]]; currentColor = stack[[1,6]];
      stack = Rest[stack];),
  (* Changing width *)
  setWidth[t_] :> (currentWidth = t)
};
pathData = Partition[Flatten[thePath],5];
widths = Thickness /@ (First /@ pathData);
lines = Line /@ (Partition[#,2]& /@ (Rest /@ pathData));
lines = Transpose[{widths, lines}] /. Thickness[t_] :> If[t<0.005,
  Unevaluated[Sequence[Thickness[t],Darker[Green]]],
  Unevaluated[Sequence[Thickness[t],Darker[Brown]]]
];

Finally, we display.

Graphics[lines]