After a long hiatus I started teaching myself Python, and as a warm-up I wrote a little lambda interpreter in it. I was blown away by how convenient and quick it was in Python, and the project has kind of snowballed to the point where I have most of my previously only vapourware “F♮” functional language up and running (albeit only under a test harness.)
The source is GPL and on GitHub: https://github.com/billhails/PyScheme
The problem: how to texture a sky sphere with procedural clouds so that the clouds appear to recede into the distance at the horizon, rather than just being wrapped around the sphere. I searched around without luck for a way to do this, before realising that it’s so easy to do, however the maths isn’t going to be obvious to everyone so I thought I’d share this quick tip.
For those who just want the node setup, here it is.
However if you’re interested here’s a brief explanation. Im assuming you’re familiar with the standard trick of using a noise texture through a colour ramp as the factor to mix between a sky texture and a clouds texture (or in this case just a plain colour). I won’t bother to explain that bit. What we want is to be able to control the scale of that cloud texture so it gets smaller as the scale value increases towards the horizon. It turns out that the value we want is just the distance from our point of view to each point on an imaginary flat horizontal cloud layer somewhere above our heads. If we take the texture co-ordinate of any point on the sky sphere as a normalised vector, then the z component of that vector is the height above ground of the end of that unit length vector. All we need to work out then is how long we would have to make that vector to touch our imaginary cloud layer.
If we say the cloud layer is 1 unit above the ground, then when z is 1, the vector is pointing straight up and the required vector length is 1. If z is 1/2, then the vector would have to be length 2 to touch the clouds, if z is 1/3, then the required length is 3 and so on. So the length of that vector is just 1/z. Plug 1/z into the scale of our noise texture and bingo, perspective clouds.
I’ve just finished reading Stephen Jay Gould’s excellent book Wonderful Life (again) and it got me thinking about random trees.
In case you haven’t read it, Wonderful Life is about the fossil bed known as the Burgess Shales which contains extrordinarily well preserved fossils of soft and hard-bodied animals from a period just after the so-called Cambrian Explosion. The Cambrian Explosion marked the period when the seas first “exploded” with an enormous range of large, multicellular animals with hard shells that preserve easily. In the 1980s a detailed re-evaluation of the fossils found in the Burgess Shales provoked a scientific revolution in paleontology, because it turns out that only a small percentage of those fossils have any direct living descendants, and many of them represent previously unknown phyla (basic types of animals.) This did not fit comfortably with the established notion of evolution as ordered progress, with the basic groups of animals established early on and forming a predictable lineage all the way from microbe to man at the pinnacle. Rather it paints the picture of extinction being the norm, and the survival of one group or another very much in the hands of chance and historical contingency. The book is not an argument against Darwinism but rather a re-evaluation of some of its finer points. Crudely put, it’s not arguing against the existance of a Tree of Life, just questioning what shape the tree is.
Anyway with that in mind, and the somewhat vague hand-drawn trees in the book leaving my curiosity piqued, I started wondering what any real evolutionary tree might look like. Of course it’s impossible to ever produce an algorithm that will accurately represent a real evolutionary sequence, so I thought to keep it very simple.
We start with a “first progenitor“. It has two choices: form two new species or die out.
Each new species has the same option at the next toss of the coin. That’s it. In perl it would look something like this:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
sub tree { my $tree = {}; my @nodes = ($tree); while (@nodes) { my $node = shift @nodes; if (rand >= 0.5) { my $l = {}; my $r = {}; $node->{l} = $l; $node->{r} = $r; push @nodes, $l, $r; } } } |
So there’s a 1/2 probability that the thing will never get started, and you’re left with a stump rather than a tree. But with 2 children, there’s only a ¼ chance that they will both die out, and if they both survive then there are 4 grandchildren, and so on. This code has a definite probability of running forever.
It turns out that if you run this a large number of times, and add up the number of each depth reached, you get a curve that asymptotically approaches zero at infinity:
The graph is normalized so the trees of depth zero come out at 0.5. The little kick at the right is those that reached the maximum depth in my test.
So what do these trees look like? I’ve given the game away by using a picture of one of them as the featured image for this post. As for generating the images, the excellent GraphViz comes to our rescue. With a little jiggery-pokery we can get the above perl code to produce a .dot file that we can feed to GraphViz and get a picture. I’ve extended the code to color nodes and vertices red if they are “survivors” (have descendents at the limiting depth) and black if they represent a species with no descendants. I’ve also changed the code to try again repeatedly until it generates a tree that reaches a limiting depth. Here’s a representative:
The limit was set at 60, so assuming 2 million years to create a species (I remember that figure from somewhere, I have a bad habit of throwing up unverified facts) this represents about 120,000,000 years of evolution from a single common ancestor. The interesting thing here I think is that the majority of branches don’t make it. Extinction is the norm, even for apparently large and flourishing branches. Apparently insignificant branches can suddenly flourish, and equally suddenly die out. I think this is close to Gould’s vision in general, if not in detail.
The other interesting thing is the huge variety of shapes. Some trees are wide, others are narrow, for example:
In this case all of the survivors share a common ancestor only a few generations (speciations) ago. This could easily be a model for the very earliest life, since the common ancestor of all current life on earth, who’s closest living relative is likely one of the Archaea, is far too complex to be a “first progenitor”.
I don’t know where I’m going with this from here, probably nowhere, but I think it’s interesting.
To finish off, here’s the full implementation of the tree generating code in case you want to try it yourself. You can pick up GraphViz from www.graphviz.org and run it from the command-line (the commands are called dot , neato , circo etc.) or via a gui.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 |
#! /usr/bin/perl use strict; use warnings; use Carp qw(confess); use constant MAX => 50; my $max = shift || MAX; sub heads { return rand > 0.5; } sub iterate { my ($sub, @nodes) = @_; while (@nodes) { my $node = shift @nodes; push @nodes, $sub->($node); } } sub mark_ancestry { my ($node) = @_; while ($node) { return if $node->{survived}; $node->{survived} = 1; $node = $node->{parent}; } } sub free { my @nodes = @_; iterate(sub { my ($node) = @_; $node->{parent} = 0; return $node->{l} ? ($node->{l}, $node->{r}) : (); }, @nodes); } sub cull { my ($tree, $depth) = @_; iterate( sub { my ($node) = @_; if ($node->{depth} == $depth) { if ($node->{l}) { free($node->{l}, $node->{r}); $node->{l} = $node->{r} = 0; } mark_ancestry($node); } elsif ($node->{l}) { return ($node->{l}, $node->{r}); } return (); }, $tree); } sub make_tree { for (;;) { my $tree = { n => 0, live => 1, depth => 1, parent => 0 }; my @nodes = ($tree); my $n = 0; my $maxdepth = 0; while (@nodes) { my $node = shift @nodes; if ($node->{depth} == $max) { cull($tree, $max); return $tree; } if (heads) { my $l = { n => ++$n, depth => $node->{depth} + 1, parent => $node, }; my $r = { n => ++$n, depth => $node->{depth} + 1, parent => $node, }; $maxdepth = $node->{depth} > $maxdepth ? $node->{depth} : $maxdepth; $node->{l} = $l; $node->{r} = $r; push @nodes, $l, $r; } else { $node->{died} = 1; } } free($tree); } } sub print_tree { my ($tree) = @_; print "digraph G {\n"; iterate(sub { my ($node) = @_; my $fillcolor = $node->{survived} ? 'red' : $node->{died} ? 'black' : 'blue'; print ' N' . $node->{n} . qq' [label="", style="filled",' . qq' fillcolor="$fillcolor" shape="circle"];\n'; if ($node->{l}) { my ($lcolor, $rcolor) = ('blue', 'blue'); $lcolor = 'red' if $node->{l}{survived}; $rcolor = 'red' if $node->{r}{survived}; print ' N' . $node->{n} . ' -> N' . $node->{l}{n} . qq' [ color="$lcolor",' . qq' arrowhead="none" ];\n'; print ' N' . $node->{n} . ' -> N' . $node->{r}{n} . qq' [ color="$rcolor",' . qq' arrowhead="none" ];\n'; return ($node->{l}, $node->{r}); } else { return (); } }, $tree); print "}\n"; } print_tree(make_tree); |
I’ve mentioned Currying and partial function application already. The idea is that given a function with more than one argument:
|
1 2 3 |
fn add (x, y) { x + y } |
if we call it with less than the arguments it expects, then it will return a function that accepts the rest:
|
1 2 3 |
def add2 = add(2,); add2(4); // 6 |
(The trailing comma is just a syntactic convention that I’ve come up with that lets the compiler know that we know what we are doing, and lets the reader know that there is Currying going on.) Now setting aside how we might type-check that, it turns out that it’s actually pretty easy to evaluate.
Normal application of a closure looks something like this (Java-ish pseudocode):
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
class Closure { private List<Symbol> fargs; private Env env; private Expr body; Expr apply(List<Expr> actual_args) { Map<Symbol, Expr> dictionary; List<Symbol> formal_args = this.fargs; while(formal_args.isNotNull() && actual_args.isNotNull()) { dictionary.put(formal_args.head(), actual_args.head()); formal_args = formal_args.tail(); actual_args = actual_args.tail(); } return this.body.eval(this.env.extend(dictionary)); } } |
For those of you that don’t know Java, List<Symbol> means List of Symbol
. And yes, we’re ignoring the possibility that we’re passed the wrong number of arguments, the type checker should deal with that.
Now if we are expecting that we might get fewer than the full set of arguments, we can instead create a new closure that expects the rest:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
class Closure { private List<Symbol> fargs; private Env env; private Expr body; Expr apply(List<Expr> actual_args) { Map<Symbol, Expr> dictionary; List<Symbol> formal_args = this.fargs; while(formal_args.isNotNull() && actual_args.isNotNull()) { dictionary.put(formal_args.head(), actual_args.head()); formal_args = formal_args.tail(); actual_args = actual_args.tail(); } if (formal_args.isNotNull()) { return new Closure(formal_args, this.body, this.env.extend(dictionary)); } else { return this.body.eval(this.env.extend(dictionary)); } } } |
Note that the dictionary that we have been building is used to extend the environment of the new closure with the values we know already, and that the formal_args we’ve been chaining down is now precisely the remaining arguments that we haven’t seen yet.
Of course this allows for silly definitions like:
|
1 2 3 4 5 |
fn add (x, y) { x + y } def anotheradd = add(,); |
But presumably our type checker (if not our grammar) would disallow that sort of thing, because there’s nothing to put a trailing comma after.
[Edit] You could alternatively add a guard clause to
apply() that says if this closure is expecting arguments and doesn’t get any, just return the original closure.
That way, something like:
|
1 |
add()()(2)()()(3) |
while still silly, would at least not be unnecessarily inefficient.
So I got the above working in F♮ easily enough, then I noticed an anomaly. The type of:
|
1 2 3 |
fn adder(x) { fn(y) { x + y } } |
is definately int → int → int, which means that the type checker is allowing it to be called like:
adder(2, 3). Why can’t I call it like that? It turns out I can:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
class Closure { private List<Symbol> fargs; private Env env; private Expr body; Expr apply(List<Expr> actual_args) { Map<Symbol, Expr> dictionary; List<Symbol> formal_args = this.fargs; while(formal_args.isNotNull() && actual_args.isNotNull()) { dictionary.put(formal_args.head(), actual_args.head()); formal_args = formal_args.tail(); actual_args = actual_args.tail(); } if (formal_args.isNotNull()) { return new Closure(formal_args, this.body, this.env.extend(dictionary)); } else if (actual_args.isNotNull()) { return this.body.eval(this.env.extend(dictionary)).apply(actual_arguments); } else { return this.body.eval(this.env.extend(dictionary)); } } } |
Assuming the type checker has done its job, then if we have any actual arguments left over then they must be destined for the function that must be the result of evaluating the body. So instead of just evaluating the body in the new env, we additionally call apply() on the result, passing in the left-over arguments.
This is pretty cool. We can have:
|
1 2 3 |
fn adder(x) { fn (y) { x + y } } |
and call it like adder(2, 3) or adder(2)(3), and we can have:
|
1 |
fn add(x, y) { x + y } |
and call it like add(2, 3) or add(2)(3).
The question arises: if we have implicit Currying, (partial function application) then do we need explicit Currying (explicitly returning a function from a function)?
The answer is a resounding yes! Consider:
|
1 2 3 4 5 6 7 8 9 |
fn bigfn (a) { if (expensive_calculation(a)) { fn (b) { cheap_op_1(b) } } else { fn (b) { cheap_op_2(b) } } } map(bigfn(x), long_list); |
We’ve only called
bigfn once, when evaluating the first argument to map, so expensive_calculation only got called once, and the explicit closure calling either
cheap_op_1 or
cheap_op_2 gets called on each element of the list.
If instead we had written:
|
1 2 3 4 5 6 7 8 9 |
fn bigfn (a, b) { if (expensive_calculation(a)) { cheap_op_1(b) } else { cheap_op_2(b) } } map(bigfn(x,), long_list); |
Then the call to expensive_calculation would get deferred until the map actually called its argument function, repeatedly, for each element of the long_list.
The Hindley-Milner Algorithm is an algorithm for determining the types of expressions. Basically it’s a formalisation of this earlier post. There’s an article on Wikipedia which is frustratingly terse and mathematical. This is my attempt to explain some of that article to myself, and to anyone else who may be interested.
The Hindley-Milner algorithm is concerned with type checking the lambda calculus, not any arbitrary programming language. However most (all?) programming language constructs can be transformed into lambda calculus. For example the lambda calculus only allows variables as formal arguments to functions, but the declaration of a temp variable:
|
1 2 3 4 |
{ var x = 10; // code that uses x } |
can be replaced by an anonymous function call with argument:
|
1 2 3 4 5 |
{ fn(x) { // code that uses x }(10); } |
Similarily the lambda calculus only treats on functions of one argument, but a function of more than one argument can be curried, etc.
We start by defining the expressions (e) we will be type-checking:
[bnf lhs=”e”]
[rhs val=”E” desc=”A primitive expression, i.e. 3.”/]
[rhs val=”s” desc=”A symbol.”/]
[rhs val=”λs.e” desc=”A function definition. s is the formal argument symbol and e is the function body (expression).”/]
[rhs val=”(e e)” desc=”The application of an expression to an expression (a function call).”/]
[/bnf]
Next we define our types (τ):
[bnf lhs=”τ”]
[rhs val=”T” desc=”A primitive type, i.e. int.”/]
[rhs val=”τ0 → τ1” desc=”A function of one argument taking type τ0 and returning type τ1“/]
[/bnf]
We need a function:
[logic_equation num=1]
[statement lhs=”f(ε, e)” rhs=”τ”/]
[/logic_equation]
where:
[logic_table]
[logic_table_row lhs=”ε” desc=”A type environment.”/]
[logic_table_row lhs=”e” desc=”An expression.”/]
[logic_table_row lhs=”τ” desc=”A type”/]
[/logic_table]
We assume we already have:
[logic_equation num=2]
[statement lhs=”f(ε, E)” rhs=”T” desc=”A set of mappings from primitive expressions to their primitive types (from 3 to int, for example.)”/]
[/logic_equation]
The following equations are logic equations. They are easy enough to read, Everything above the line are assumptions. The statement below the line should follow if the assumptions are true.
Our second assumption is:
[logic_equation num=3]
[assumption lhs=”(s, τ)” op=”∈” rhs=”ε” desc=”If (s, τ) is in ε (i.e. if ε has a mapping from s to τ)”/]
[conclusion lhs=”f(ε, s)” rhs=”τ” desc=”Then in the context of ε, s is a τ”/]
[/logic_equation]
Informally symbols are looked up in the type environment
.
[logic_equation num=4]
[assumption lhs=”f(ε, g)” rhs=”τ1 → τ” desc=”If g is a function mapping a τ1 to a τ”/]
[assumption lhs=”f(ε, e)” rhs=”τ1” desc=”and e is a τ1“/]
[conclusion lhs=”f(ε, (g e))” rhs=”τ” desc=”Then the application of g to e is a τ”/]
[/logic_equation]
That is just common sense.
[logic_equation num=5]
[assumption lhs=”ε1” rhs=”ε ∪ (s, τ)” desc=”If ε1 is ε extended by (s, τ), e.g. if s is a τ”/]
[conclusion lhs=”f(ε, λs.e)” rhs=”τ → f(ε1, e)” desc=”Then the output type of a function with argument s of type τ, and body e, is the type of the body e in the context of ε1“/]
[/logic_equation]
This is just a bit tricky. We don’t necessarily know the value of τ when evaluating this expression, but that’s what logic variables are for.
new which returns a fresh type variable each time it is called.eq which unifies two equations.We modify our function, part [4] (function application) as follows:
[logic_equation num=6]
[assumption lhs=”τ0” rhs=”f(ε, e0)” desc=”If τ0 is the type of e0“/]
[assumption lhs=”τ1” rhs=”f(ε, e1)” desc=”and τ1 is the type of e1“/]
[assumption lhs=”τ” rhs=”new” desc=”and τ is a fresh type variable”/]
[conclusion lhs=”f(ε, (e0 e1))” rhs=”eq(τ0, τ1 → τ); τ” desc=”Then after unifying τ0 with τ1 → τ, the type of (e0 e1) is τ.”/]
[/logic_equation]
That deserves a bit of discussion. We know e0 is a function, so it must have a type τa → τb for some types τa and τb. We calculate τ0 as the provisional type of e0 and τ1 as the type of e1, then create a new type variable τ to hold the type of (e0 e1).
Suppose e0 is the function length (the length of a list of some unspecified type τ2), then τ0 should come out as [τ2] → int (using [x] to mean list of x
.)
Suppose further that τ1 is char.
We therefore unify:
[logic_equation]
[statement lhs=”{{{τ2}}}” op=”→” rhs=”int“/]
[statement lhs=”{{{char}}}” op=”→” rhs=”τ”/]
[/logic_equation]
Which correctly infers that the type of (length [‘c’]) is int. Unfortunately, in doing so, we permanently unify τ2 with char, forcing length to have permanent type [char] → int so this algorithm does not cope with polymorphic functions such as length.
The full score and separate parts as PDF are free for download. If you make anything of them I’d be pleased to know:
