representations, permutations, visualisations

One of the things I’m interested in is evolutionary algorithms (EAs), and how to make them better.  An EA takes a population of “genomes”, “mutates” (changes a little) and selects (based on “fitness”), and mutates and selects, and ... until a suitably fit answer is found.

A recent advance has been the introduction of “evo-devo” algorithms.  (I’m putting all this biological terminology in scare quotes, because by the time the relevant process has been translated into a computer algorithm, it is so far removed from its biological inspiration as to make a biologist wince, or even exclaim in outrage.)  Evo-devo puts a distance between the genome (the representation that gets mutated) and the “phenotype” (the representation that gets selected based on fitness).  This can help the algorithm’s performance, by allowing simple easily mutatable genomes develop into complex structured phenotypes.

A colleague of mine at York, Jillian Miller, is the inventor of such an algorithm, Cartesian Genetic Programming (CGP).  The genome is a string of numbers (numbers are easy to mutate a little bit).  The string is then interpreted as a network phenotype.  The network itself has inputs and outputs, so is a form of program.

Now, let’s consider permutations.  A permutation of the numbers 1 to N is these numbers in some specific order.  So the permutations of 1 to 3 are: (1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), (3,2,1).  A string of length N has N! (N factorial) permutations.  N! grows very fast; while 5! = 120, 10! = 3,628,800, and 100! > 10¹⁵⁷.

Permutations are common in computer science.  One classic use is in the Travelling Salesman Problem: given a bunch of N cities, find the shortest path through all of them.  That is, find the permutation of 1 to N that gives the shortest path.  Given there are N! such permutations, clearly we don’t want to try them all.  Although exact algorithms that are essentially more efficient than trying all possibilities aren’t known (and it is strongly suspected that there aren’t any), there are algorithms that come up with very good approximate (nearly shortest path) answers most of the time.  EAs are one such class of algorithms: breed for fitter (shorter) paths.

Permutations as genomes are a bit tricky, though.  A permutation has structure: it must contain all the numbers from 1 to N, and each only once.  So you can’t mutate a single entry: you have to swap two entries, or do something else that maintains the permutation structure. “Crossover” is even harder: how do you take half of one permutation, half of another, and combine them into a valid permutation?  There are various techniques, but they are not very pretty.

Using an evo-devo approach to generate a permutation seems even harder: how do you ensure that your developed system is a valid permutation?  So, for example, with CGP we can have a list of outputs, but how do we ensure that this list is a valid permutation?  (Having a single output that is already a permutation merely moves the problem back inside the network somewhere.)

We need a further representation and development step that is guaranteed to produce a permutation.  Rather than try to get the network to produce a permutation immediately, let’s break it down into two steps: the network produces a list of numbers, then that list has to go through a further interpretation step to form a permutation.  Julian came up with an idea of how to do this: given a list of (say) real numbers (easy to produce with CGP), just sort them into ascending order.  The correspondingly sorted list of the indexes gives the required permutation.  Voilà!

A string of numbers is interpreted as a network, which outputs a real vector, which when sorted yields a permutation

What is happening here is easy to visualise using a technique called parallel coordinates.  A list of N real numbers can be thought of as a vector in N-D space.  But N-D space is hard to visualise if N > 3.  (I find it pretty hard to visualise even when N = 3.)

It’s hard to visualise a lot of dimensions this way

Parallel coordinates do what it says in the name: instead of drawing the N dimensions orthogonal to each other (rapidly running out of ways to do this in our 3D physical space), draw them parallel to each other.  It’s easy to draw lots of parallel lines.  Now plot the N-D point (x₁,x₂,...x_N) as follows: plot the point x₁ on axis 1, the point x₂ on axis 2, and so on, then joint these points together with a line.  The line in the parallel coordinate plot represents the point in N-D space.

parallel coordinates view of a single N-D point

We can use these parallel coordinated to visualise how a vector of real numbers can represent a permutation by its components being sorted into ascending order.

(top) a vector of 20 real numbers, a 20-D point, drawn in parallel coordinates; (bottom) the same vector, with the parallel axes ordered so that the components are in increasing order: the sorted axis indexes are the permutation represented by the N-D point. 

The Python/numpy code that generated these plots is:

N = 20
P = range(N) # indexes
V = rand(N)  # random vector

# plot unsorted vector
for dim in range(N):
    ax1.plot([dim, dim], [0, 1], '0.5', linewidth=0.25)
    ax1.text(dim, -0.2, str(P[dim]), ha='center', fontsize=32)
ax1.plot(range(N), V, '.k', markersize=20)
ax1.plot(range(N), V, 'k')

# sort, and plot sorted vector
P = argsort(V) # sort indexes
V = sort(V)    # sort vector (for plotting)
for dim in range(N):
    ax2.plot([dim, dim], [0, 1], '0.5', linewidth=0.25)
    ax2.text(dim, -0.2, str(P[dim]), ha='center', fontsize=32)
ax2.plot(range(N), V, '.k', markersize=20)
ax2.plot(range(N), V, 'k')

Note that the code that generates the permutation is the single line P = argsort(V): the rest is just plotting code.

Here I started from a random vector, rather than the non-random output of some CGP network.  Sorting a random vector is one way to construct a random permutation, but as far as Julian and I can tell from the literature, this CGP use for representing evolved, non-random permutations isn’t standard.  Julian has been using it for several years in his module on evolutionary algorithms, and will be publishing a paper on some results next year.

Read More

goldfish brain

For the last eleven years and more, I’ve been driving between home near Cambridge, and work in York, nearly every week. I know the 350 mile round trip rather well by now. Despite this, I’ve taken to having the SatNav on during the journey. That’s because there can often be delays due to accidents, or road closures due to roadworks, and a little advance warning can let me choose a route more useful than the suggested diversion.

Last week, the inevitable happened. I was tootling up the A1 to York, when the SatNav advised me: “Road closed at J36 due to roadworks. Calculate alternative route?” Yes.

Whirr, whirr, whirr. “Confirm new route?” The map is very small, but I could see this new route at least approximated the M18 diversion I had already intended to take. So, yes, let’s take that one.

Pause. “Road closed at <somewhere I’d never heard of> due to roadworks. Calculate alternative route?” Okay, that’s a nuisance, so yes, calculate new route.

Whirr, whirr, whirr. “Confirm new route?” The map is very small, but I could see this new route looked close to the original A1 route.  But, hey.  Yes, confirm new route.

Pause. “Road closed at J36 due to roadworks. Calculate alternative route?”  Aargh! it had put me back on the original closed route.  I could guess what would happen next, but I tried it anyhow.  Yes, calculate alternative route.

Whirr, whirr, whirr. “Confirm new route?” Back on the M18 diversion. Yes, confirm. Pause. “Road closed at <somewhere I’d never heard of until a few minutes ago> due to roadworks. Calculate alternative route?”

Clearly the system has a diversion buffer of size one.  I decided to take the M18 diversion. Once I was on that bit of road, the system the quite happily found me a diversion around the diversion.

(left) intended route, but the A1 was closed around Doncaster; (mid) SatNav suggested alternative, but then found to be closed somewhere south of Selby; (right) final SatNav suggested route, once commited to the M18 east of Doncaster

very #FirstWorldProblems

Read More

normality has been restored...

...partially, anyway.

When my home machine died, and I upgraded to Windows 7, my incredibly ancient 16-bit Smalltalk/V would no longer work.  As this was what I used to generate the book review pages on my website, it was a bit of a disaster.

So over the past several weekends, I've been re-implementing my book list program, from scratch, in Python, and the new pages have finally gone live.  (There are still broken links from the non-automatically generated pages; I'll fix those next.)

This was my first major foray into Python.  Three years ago, I wrote my first ever Python program -- a Game of Life implementation, using Gosper's hashlife algorithm -- but then didn't touch it again until recently.  Instead, I was using Matlab for all my number-crunching needs.  However, I needed something a bit more friendly for teaching some non-computer science students a little bit about simulation.  So before Christmas I wrote a small flocking program in Python, and a few little bits and bobs.  Then I learned about NumPy, which does a lot of the matrix stuff that Matlab does, so I thought I'd give it a try, by reimplementing some Random Boolean Network code I'd originally written in Matlab.

I was really impressed, so decided to do some more NumPy programming, to learn it better.  I wrote a few programs to support some practicals for a module I teach: a quantum random walk simulator, an iterated function system, and a 1D cellular automaton.  And after that exposure, I decided I was deeply in love with Python.

But the book list program is my biggest yet, by a long way.  All the previous programs had a core of essentially about 10 lines of code (each of those lines doing a lot, mind you!), each written in an hour or so.  The book list program is several hundreds of lines (many of which also do rather a lot), split over 5 files and several classes.  I also decided to take the opportunity to reorganise and tidy up the structure of my web pages, which involved writing a few utilities, to rename a bunch of files, to write a new table in the database, that sort of thing.  File handling. Database access. All amazingly simple in the marvellous language that is Python.

Although, frankly, I couldn't have done it without the wonderful Stackoverflow forum.  The Python documentation is useful if you know how to do something, and just need to check the syntax or method name.  It's not so helpful if you don't know how to do something. And Python has its own way of doing things.  But every time I typed into Google: "python how to  XXX", for a wide range of XXX, there was Stackoverflow with the answer -- often several answers!  Thanks, guys.

When I started writing this, I was consulting online resources for nearly every line.  But now, I'm happily typing things like

reviews = [bk for bk in have_reviews if bk.isNFBook()][0:10]

without looking things up any more.

And I'm now looking round for my next Python project...

Read More

faceless

I was just uploading my photo to my Google profile.  The drag-and-drop interface is nice (although it's very irritating to be told of the minimum 250x250 pixel size only after uploading a smaller pic).  However, I was a little disconcerted to see:

No, the problem isn't that I have acquired a large red halo!  That's just there to highlight the surprising warning box:

I'm not sure whether to be concerned at my lack of face, or to be pleased that I'm immune to Google's face recognition software...

Read More

RBNs with NumPy, sorted

I've been using Python for a little while now, and love the ease of programming in it.  I also use Matlab, which is wonderful for programming scientific things, particularly with arrays.  But Matlab is expensive, so I only have access to it at work.  Python is free.

I heard that NumPy, the numerical package for Python, had Matlab-like array operations, so thought I'd give it a try.  This weekend I finally had some time (that is, I needed a displacement activity from marking), so I gave it a go.  I decided to do a comparison of something I'd already implemented in Matlab: a Random Boolean Network (RBN) tool.

RBNs were invented by Stuart Kauffman as a simplified model of gene regulatory networks. They have fascintating properties for so simple a construction.  An RBN has the following components:

• +++N+++ binary-state nodes: +++n_1 .. n_N+++
• each node can be off or on (in state 0 or 1), +++s_1 .. s_N+++
• each node has +++K+++ input connections from +++K+++ randomly chosen different nodes in the network, +++c_{11} .. c_{1K} .. c_{N1} .. c_{NK}+++
• each node has a randomly chosen boolean function of +++K+++ variables, +++b_1 .. b_N+++

An example of an +++N=5, K=2+++ RBN is

Here the node colours represent the different boolean functions +++b_i+++, and the numbers label the nodes from +++0 .. N-1+++.

You start the network in some initial state of the binary nodes.  Each timestep each node receives the state of its +++K+++ connected neighbours, combines them with its boolean function, and sets its next state to that value:
$$s_i(t+1) = b_i(s_{c_{i1}}, s_{c_{i2}}, ... s_{c_{iK}})$$
The marvellous thing about +++K=2+++ RBNs is that, despite being set up to be as random as possible, and having a total of +++2^N+++ possible states they could be in, they rapidly settle down into an attractor cycle of length +++O(\sqrt N)+++.

This establishment of order from seeming randomness is fascinating, but really needs to be demonstrated to be appreciated.  Hence NumPy.

Here's the code:

import matplotlib.pyplot as plt
from numpy import *

K = 2       # number of connections
N = 500     # number of nodes, indexed 0 .. N-1
T = 200     # timesteps

Pow = 2**arange(K) # [ 1 2 4 ... ], for converting inputs to numerical value
Con = apply_along_axis(random.permutation, 1, tile(range(N), (N,1) ))[:, 0:K]
Bool = random.randint(0, 2, (N, 2**K))

State = zeros((T+1,N),dtype=int)
State[0] = random.randint(0, 2, N)
for t in range(T):  # 0 .. T-1
    State[t+1] = Bool[:, sum(Pow * State[t,Con],1)].diagonal()

plt.imshow(State, cmap='Greys', interpolation='None')
plt.show()

There's essentially only three lines doing much substantive, which is the joy of working directly with arrays: no fiddly, wordy iterations.  Con holds the random connections, Bool holds the random functions, and the loop over t calculates the next State each timestep.

I dare say there's more elegant ways to do this, but I am still learning NumPy's capabilities. But what exactly is going on here?

For the Con array we need to choose K random inputs for each node.  These need to be distinct inputs, so we can't just choose them at random, because there might be collisions.  We could keep choosing, and keep throwing away collisions, but there's another way to do it:

• range(N) gives a list [0, 1, .., N-1].  Let's take N=5 here
• tile(...) makes an array of 5 stacked copies of this:
• [ [ 0 1 2 3 4 5 ]
  [ 0 1 2 3 4 5 ]
  [ 0 1 2 3 4 5 ]
  [ 0 1 2 3 4 5 ]
  [ 0 1 2 3 4 5 ] ]
• apply_along_axis(random.permutation, ...) applies a random permutation to each row individually
• [ [ 2 0 1 3 5 4 ]
  [ 0 2 3 4 5 1 ]
  [ 1 3 2 5 0 4 ]
  [ 3 1 4 2 0 5 ]
  [ 4 2 0 5 1 3 ] ]
• ...[:, 0:K] takes the first K items from each row.  Here K = 2
• [ [ 2 0 ]
  [ 0 2 ]
  [ 1 3 ]
  [ 3 1 ]
  [ 4 2 ] ]

This gives the array of node connections: node 0 has inputs from node 2 and itself; node 1 has inputs from node 0 and node 2, and so on.  See the figure earlier.

Bool has the +++N+++ random boolean functions.  Here each function is stored as a lookup table: a list of +++2^K+++ ones and zeros.

For the State update

• State[t,Con] gets the inputs from the connections
• sum(...) converts this array of ones and zeros into an index +++0..2^K-1+++
• Bool[: ...].diagonal() looks up the next state value from this index 

And that's it!  We can use this to plot the time evolution of an RBN.

The string of nodes is drawn as a horizontal line at each timestep, and  time increases down the page. You can see it has some random structure for the first few timesteps, then rapidly settles down into regular behaviour.

Well, it's not that easy to see the regular behaviour.  It's a bit of a jumble really.  We can do better.

The problem is, since an RBN is random, there's appears to be no obvious order to write down the nodes.  The picture above uses the order as first given, which is ... random.  However, an RBN does have structure; it has a "frozen core" of nodes that settle down into a "frozen" state of always on, or always off.  If we sort the nodes by their overall activity, it highlights the structure better.

So, there's a little bit of extra code, sitting just before the loop updating the state.

SRun = 5     # sorting runs
ST = 200     # sorting timesteps
State = zeros((ST+1,N),dtype=int)
Totals = State[0]

for r in range(SRun):
    for t in range(ST):
        State[t+1] = Bool[:, sum(Pow * State[t,Con],1)].diagonal()
        Totals = Totals + State[t+1]
    State[0] = random.randint(0, 2, N) # new initial random state
    
Index = argsort(Totals)    # permutation indexes for sorted order
Bool = Bool[Index]         # permute the boolean functions
Con = Con[Index]           # permute the connections

InvIndex = argsort(Index)  # inverse permutation
Con = InvIndex[Con]        # relabel the connections

This extra code runs the RBN several times (from different initial conditions, each potentially leading to different attractor cycles involving different patterns of node activity), totalling up the number of times each node is active.  Sorting this array puts more inactive nodes towards the start, and more active nodes towards the end.  argsort() doesn't return the sorted array, however; it returns a permutation of the indexes corresponding to this sort.  This Index array can then be used to sort the Con and Bool arrays into the same order.  This results in something like:

Having done this, we need to relabel the nodes and the connection indexes.  This requires using the inverse sort permutation, InvIndex.  So we get

Running the modified code gives a much clearer picture of the RBN's dynamic behaviour:

So, after all this, what do I think of NumPy?

It's excellent.  Everything I needed (random permutations, sorting, applying functions across arrays, indexing arrays with other arrays, whatever), it's all there.  The code produced is very compact. (So compact that I've commented it quite liberally in the source file.)

The online documentation is mostly adequate, and whenever I puzzled over how to do something, a quick Google usually got me to a forum where my question had already been answered.

NumPy also has a big advantage over Matlab (in addition to the price!).  With Matlab, many functions and operations (such as array indexing) can be applied only to array literals, not to array expressions.  This makes it hard to build up compound operations without having to have a lot of intermediate variables.  With NumPy, you can just build up the expression in one go.  That makes for a more natural style of programming (although I suspect it could also make for some spectacular write-only code).

Anyhow, I'm please with my experiment, and will be delving further into NumPy in the future.

Read More

ODE to a Petri net

Writing ordinary differential equations (ODEs) to model various natural world processes comes more readily to some than to others.  And once written, it can take some effort to pick apart their real-world meaning.  It would be nice to have a more visual form.  Petri nets are one such approach.

Epidemics

For example, consider the simple SIR model of epidemic infection.

• +++S+++ : those uninfected, but susceptible -- their number is reduced as they become infected, at a rate proportional to the number of susceptibles and the number of infected
• +++I+++ : those infected -- their number is increased by susceptibles who become infected;  it is also reduced as those infected recover, at a rate proportional to the number of infected
• +++R+++ : those recovered -- their number is increased by those who were infected recovering

For those happy with ODEs, it is straightforward to write down a set of coupled equations to model this:

$$\dot{S} = - i S I $$

$$\dot{I} = i S I - r I$$

$$\dot{R} = r I$$
In addition to natural language (the bullet list) and maths (the ODEs), there is another language useful to explain and understand models: diagrams.

For example, we could draw a simple state transition diagram to show the movement from susceptible to infected to recovered:

This captures some of the information, but not all of it.  A (continuous) Petri net can do better:

The circles are called "places", and represent the "things" involved: here the susceptibles, the infected, and the recovered.  The rectangles (the colours aren't significant) are called "transitions", and represent how the things in "input" places get transformed into things in "output" places.

• transition +++r+++: an infected comes in, and a recovered comes out
• transition +++i+++: a susceptible and an infected come in; two infected come out (the two infected outputs are the newly infected, and the original infecter)

This diagram has enough information in it to reproduce the original ODEs.  We have one ODE per place, with the terms being given by the transitions feeding that place.

• place +++R+++. +++R+++ has only one transition feeding it: transition +++r+++.  It is feeding into +++R+++ at a rate proportional to all the inputs to +++r+++: here just +++I+++.  Hence +++\dot{R} \propto I+++.  If we call the constant of proportionality (the rate constant) +++r+++, we get +++\dot{R} = r I+++
• place +++S+++. +++S+++ has only one transition, +++i+++, removing stuff from +++S+++.  It is removing at a rate proportional to all the inputs to +++i+++: here +++S+++ and +++I+++.  Calling the rate constant +++i+++, we get  +++\dot{S} = - i S I+++ (the minus sign is there because we a removing from +++S+++, and it is conventional to keep the rate constants positive)
• place +++I+++.  +++I+++ has two transitions: +++i+++ both feeding it (two input arrows) and removing from it (one output arrow), and +++r+++ removing from it.  We get +++\dot{I} = i S I - r I+++

Thus we have recovered the original equations.

We can write this focussing on the transitions, in a more algorithmic way (algorithms, or pseudo-code, are yet a further language we can use to explain, describe, and define processes).

• initialise the rates of change to each place to be +++0+++.  +++\dot{S}, \dot{I}, \dot{R} := 0+++
• transition +++r+++.  This is removing stuff from place +++I+++ and adding it to +++R+++ at a rate +++rI+++. Update the output place +++R+++ and input place +++I+++ appropriately: +++\dot{I} {-}{=}\ rI+++, +++\dot{R} {+}{=}\ rI+++
• transition +++i+++.  This is removing stuff from place +++I+++ and from place +++S+++ at a rate +++iSI+++ and adding it to +++I+++ at a rate +++2iSI+++ (from the two input arrows).  Hence there is a net input to place +++I+++ at a rate +++iSI+++.  Update the output place +++S+++ and net input place +++I+++ appropriately: +++\dot{I} {+}{=}\ iSI+++, +++\dot{S} {-}{=}\ iSI+++

We can write this as a general algorithm:

  for each place Pi
    Pi_dot := 0
  for each transition Ti
    let Pin = < Pin_1, ... , Pin_n > = list of n places,
                                       one for each input arrow of Ti;
        Pout = list of m places, one for each output arrow of Ti;
        t = Ti x Pin_1 x ... x Pin_n
    for each Pi in Pin
      Pi_dot -= t
    for each Pi in Pout
      Pi_dot += t

So now we have a diagrammatic form, and an ODE form, that are equivalent, and an algorithm to translate on to the other.  This is useful, because we can use them interchangeably, without risk of losing information.  In particular, notice how explanation accompanying the Petri net focusses on what is happening in the transitions, whilst that for the ODE form focusses on what is happening to the places.  Having different forms of explanation can be useful in different circumstances (modelling, communication, modification, validation, calculation, etc).

Catalysis

Although that all works well, the handling of the infecter in the +++i+++ transition seems a bit unnatural: infecter goes in, infecter comes out, resulting in an addition and subtraction of this rate.  The infecter is needed for the transition, and affects the rate of the transition, but is not itself changed by the transition.  In chemistry, this is called a catalyst, and there is some special Petri net syntax for it.  We can draw the SIR Petri net above equivalently as:

Here the dashed arrow means that +++I+++ is a catalyst: it is needed for the transition, but is not consumed by the transition. Hence there is now only one arrow out to +++I+++: since the catalyst wasn't consumed, it doesn't need to be replaced; the remaining single arrow represents the newly infected.

The algorithm needs a little  modification:

  for each place Pi
    Pi_dot := 0
  for each transition Ti
    let Pin = < Pin_1, ... , Pin_n > = list of n places, 
                                       one for each input arrow of Ti;
        Pcat (sublist of Pin) = list of catalytic input arrow of Ti;
        Pout = list of m places, one for each output arrow of Ti
    let t = Ti x Pin_1 x ... x Pin_n
    for each Pi in Pin \ Pcat
      Pi_dot -= t
    for each Pi in Pout
      Pi_dot += t

So the catalytic arrows still contribute to the functional form of the overall rate +++t+++, but not to the changes to the specific places.

Logistic equation

Possibly the simplest bounded growth model in biology is the logistic equation:
$$ \dot{N} = rN(1-N/K)$$where +++r+++ is the growth rate, and +++K+++ is the carrying capacity (so when +++K=N+++, +++\dot{N}=0+++).

This can be drawn as an equivalent Petri net:

• transition +++r+++ (birth): one in, two out
• transition +++rK+++ (competition): two in, one out

These two transitions can also be shown in a simpler catalytic form (if maybe not with the same intuition as before):

• transition +++r+++ (birth): one "catalyses" the birth of the other
• transition +++rK+++ (competition): one "catalyses" the death of the other

Lotka-Volterra predator-prey

The simple predator-prey model, usually cast as rabbits and foxes, has rabbits being born, predated on by foxes to produce more foxes, who then die.  A simplistic version of this might be:

• transition +++b+++ (birth): one rabbit in, two rabbits out
• transition +++d+++ (death): one fox dies
• transition +++p+++ (predation): one fox and one rabbit in, two foxes out

This however has a problem: it has a new fox produced every time a rabbit is eaten.  Real foxes need more food than this to reproduce.  We can't solve the problem by changing the rate +++p+++, as this affects the consumption of rabbits and production of foxes equally.  What we really need is for the consumption of a rabbit to produce a bit of a fox.  We can do this by adding a separate rate to the arrow:

• transition +++p+++ (predation): one fox and one rabbit in, one plus +++\epsilon+++ foxes out

The algorithm needs to be updated to multiply the rate by the weight of the arrow before adding/subtracting  as appropriate.  This then yields the familiar equations:

$$\dot{R} = R(b-pF)$$

$$\dot{F} = F(-d+p\epsilon R)$$

If we use the catalytic form, the diagram simplifies to:

• transition +++p+++ (predation): one fox "catalyses" the transformation of a rabbit into +++\epsilon+++ of a new fox

Lotka-Volterra competition

The simple competition model, usually cast as rabbits and sheep, has rabbits and sheep being born and dying following their own logistic equation, and also competing with each other for resources.

• transition +++rb,rs+++ (birth): one +++x+++ "catalyses" the birth of the next +++x+++
• transition +++rc,sc+++ (competition): one +++x+++ "catalyses" the death of another +++x+++
• transition +++rs+++ (rabbit/sheep competition): one sheep and one rabbit in, a proportion of each out

$$\dot{R} = R(rb-rcS) + (pr-1)rs RS$$

$$\dot{S} = S(sb-scS) + (ps-1)rs RS$$

Combining predator-prey and competition

The competition example has two logistic "subnets", showing how these diagrammatic forms can be readily combined.  So, for example, we could easily add some foxes to the brew:

$$\dot{S} = S(sb-scS) + (ps-1)rs RS$$

$$\dot{R} = R(rb-rcS) + (pr-1)rs RS -pRF$$

$$\dot{F} = F(-d+p\epsilon R)$$

If the foxes also worried the sheep, the diagram would get messier, but would still visually represent the relationships between the different components.

Diagrams v cartoons

Pictures can be very helpful at getting across ideas, but they have their problems if they are ambiguous, incomplete, or otherwise open to misinterpretation.

The (continuous) Petri nets shown here have a formal meaning: they can be translated into equivalent ODEs. They are not informal "cartoons", merely sketching some part of the meaning.  There is an algorithm from diagrams to equations, and it is possible to build a tools allowing the manipulation of diagrams and equations as two different "concrete syntaxes" of the same underlying model.  Which syntax to use depends on what you are doing: it truly is the best of both worlds.

Acknowledgments

• I first came across the formal link between continuous Petri nets and ODEs on Alexi Sharov's web site
• I was reminded of using Petri nets to model population dynamics on reading David Tanzer's guest post on the Azimuth site
• I drew the diagrams in graphviz
• The maths is formatted with LaTeX and displayed using MathJax

Read More

in praise of waste paper

I have a stack of "rough paper": paper printed on one side but no longer needed, that I keep by my desk for scribbled calculations and diagrams, for when I am thinking.  I've never really considered this, except occasionally to think that I should actually scribble these notes into my daybook as I go, so that they are part of my technical log.  That has always felt wrong to me for some reason.  Now I know why: Venkatesh Rao has the same experience, but understands it:

I can’t really work with ruled paper or blank two-sided paper. Both are pristine resources that I feel guilty about wasting with my manic scribbling of mostly useless thoughts. ...
... the reason I work best with one-sided printed paper is that it has already been bad-wasted. It’s going to the recycling bin anyway. So it is a particularly liberating medium to work with. More than liberating. You can feel virtuous because you are effectively redeeming bad waste.

Yes, exactly this.  I can draw a small diagram in the middle of an A4 sheet, dislike it, discard it, and draw again, and again, with no feelings of guilt (like Venkat, I use other people's discarded paper).  This freedom to redo things is very liberating.  If I was drawing in my daybook, or using some pristine pad of blank paper, I would feel constrained to be more frugal in my explorations.  As it is, I give the paper no thought whatsoever (except maybe a small glow of productivity when I discard a great wodge of it in the recycling bin at the end of the day: look at all that work I must have done!).

Venkat's post is about waste-enabled creativity: "civilizations are defined by the resources they can waste".  He talks about exploration-exploitation tradeoffs, where being able to waste resources allows cheaper, and hence more, trial-and-error exploration, and hence the discovery of potentially better solutions to exploit.

It's not just civilisations.  Biological life can be similarly defined.  Computer Scientists often (incorrectly) assume that biological evolution is a great optimisation algorithm.   This is incorrect for two reasons: (1) it's not "great"; it's wasteful: millions are born to give a varied pool from which a few fit are selected; (2) it's not an optimiser.  Biological life "wastes" enormous quantities of resource (a million cod are spawned so that two may survive to reproduce); it explores the vast evolutionary landscape like mad, and is fantastically creative as a result.

We Computer Scientists love bio-inspired algorithms, because of the (assumed) robustness of biological organisms.  But if only, we often wish, especially when we are being "good" frugal engineers, if only we could have all those lovely biological properties without all that waste.  Hmm.  Maybe in our attempt at software engineering and Artificial Life, the problem is that we are not wasting enough, and that we are wasting the wrong things. 

Read More

rough fonts

I was reading the Boing Boing post about fonts that look like graphs. More interesting was the link to "postscript hacks to render slightly differently every time". I followed up, and discovered the Beowolf font family, with parameterised roughness, caused by adding some randomness to the PostScript rendering process:

Beowolf family, from MyFonts 

What's great about this site is the number of different kinds of text you can choose to display the various fonts.  I find the mid-rough R22 style most pleasing to my eye. 

A 2008 post on the Typographica blog gives some interesting historical background to this family.  (There's a blog all about fonts!?!  I may never get anything useful done again...)

This reminds me of "roughness", which is one of Christopher Alexander's 15 "generative properties", from his Nature of Order work, used to produce architectural structures that exhibit "wholeness".  We did some work a while back trying to formalise a few of the other properties, through BlobWorld, and were interested to discover an optical illusion that led to an appearance of roughness emerging.  Rough brickwork, or rough wall alignment, okay.  But rough fonts?  Neat.

Read More

a stopped clock is never right

A quarter past six? Or is it?

There's a saying that

Even a stopped clock is right twice a day.

This is used to mean: even something completely unreliable can (accidentally) sometimes be right.  The saying is sometimes cast as a paradox

A stopped clock is better than a clock an hour slow, because a stopped clock is right twice a day, yet a clock an hour slow is never right!

I want to explain how, in fact, a stopped clock is never right, and a clock an hour slow is always right.  This requires us to think of a clock as a simple computer, computing the current time, and ask ourselves, how can we tell the current time from the output of its computation?

A computation has three steps

1. initialisation: set up the computer to perform the task of interest
2. operation: the computer does its thing
3. finalisation: read off the answer from the computer

(Don't blame me for the step names; I didn't invent them!  For those of you who are interested, this terminology comes from computational refinement theory.)  For a clock, these steps are instantiated as:

1. initialisation: set the clock to the current time
2. operation: the clock does its thing, marking off the passing moments
3. finalisation: read off the (now later) current time from the clock

It's nearly ten past ten. Or is it?

Notice how the finalisation step is non-trivial.  The clock doesn't output "the time": it displays an output that requires some effort to be interpreted as the time.  To read the time from my analogue wristwatch (yes, I still use a wristwatch, and yes, it has an analogue display), I have to convert the positions of the hands, relative to a standard vertical (12 o'clock!) position, into a time.  This takes a (small amount of) skill: I can remember being taught how to "tell the time", that is, read an analogue clock face, by my aunt when I was about five. Even to read the time from a digital face requires some processing: to convert the displayed pattern of LED segments, or of pixels, into characters (ie, to read the  displayed pattern as characters), and interpret those characters ("12:30", say) as a time ("half past twelve").

This finalisation step is not the only one that can be applied, however.  This is the key step in the argument.  I realised this when I was attending a conference in Toulouse in 1999, and my watch was "broken".  It hadn't stopped, but I couldn't change the time (I couldn't reinitialise it), so it was an hour slow (stuck on UK time). That is, when I interpreted its output using the conventional finalisation, the time I got was off by one hour.  Given we started with "a clock an hour slow is never right", I could nevertheless use my watch to tell the correct time.  How?  (The answer will be obvious to anyone who has used a sundial during daylight savings time.)  By applying a different finalisation, one appropriate to its actual initialisation.  Here's the setup:

1. initialisation: set my watch to the current UK time, so an hour behind the current French time
   
2. operation: the watch does its thing, ticking off the passing moments
   
3. finalisation: read off the (now later) time displayed by my watch, and add one hour

Voila! My watch was computing the correct French time, provided I finalised it correctly, that is, provided I correctly interpreted its output.  Well, you might say, but how did you know to add the hour?  Because I was the one who initialised it: I was the one who set up the computation. Other people looking at my watch would be confused, because they would be applying a different finalisation: the conventional one.  But it is merely a convention (established to make it convenient to use clocks other than ones that you have set yourself). In truth, you cannot tell the time looking at a clock unless you have some extra information: what finalisation you need to apply to interpret its display as a time. In practice, applying the conventional finalisation works, most of the time.

By using even more powerful finalisations, we can compute the time using even more faulty watches. For example, if I have a watch that loses a minute every hour, I can still use it, by adding the correct number of minutes back on when I read the time.  It is the combination of operation and finalisation that gives the resulting computation.

So how about the stopped clock?  Can you use it by applying an even more powerful finalisation? No. There is no finalisation that allows you to read off the correct time. The clock is performing no operation, it is not marking the passing time, so in order to get the desired computation from it, all the work would have to be done in finalisation alone, which would require using another clock!  The stopped clock is never right, because there is no finalisation: no way to interpret the display.  Your fortuitously looking at it when its display shows the current time does not make it right, not even coincidentally, because you have no way of interpreting its output.

However, there is something that a stopped clock has computed: the time at which it stopped (subject to applying the correct finalisation to its display, the one that would have been used when it was working).  This computation is the staple of many a TV cop show to tell the time of death of the newly discovered corpse with a conveniently smashed wristwatch (and the conventional finalisation being the wrong one is a cunning red herring in several detective novels).

Read More

Regress bars

I was copying a large set of files, when I noticed some peculiar behaviour of the "progress" bar.

I'm used to the Minutes Remaining indicator usually being pretty meaningless (I've seen it go from days, to hours, to minutes, to seconds, over the course of seconds, and it often leaps up by an order of magnitude if I start doing something else on the computer). But this time it seemed even worse than useless.  I watched it for a few moments, then felt compelled to grab a pen and paper.  (I used an antiquated technology to record the data, as I didn't want to interfere with the copy process. These fluctuations were all its own work.)

The whole exercise covered a period of a couple of minutes (that is, the line representing the actual time remaining would be almost indistinguishable from the x-axis on this plot).  The final estimate (of 15 seconds remaining, up from the previous 10 seconds) stayed there for about three seconds, and the copy was finally finished.

Well, as they say, "Prediction is difficult, especially about the future".  But one might hope that the algorithm used gave better predictions as time progressed...

Read More

name that tune

I have great difficulty identifying tunes, even if I know them well. I have no way of asking for help, since all notes I sing are the same, and flat, even though I can hear the tune perfectly well inside my head.

So, many years ago I was excited to hear of a book that would help: it lists all tunes in order of whether their notes go up, down, or stay the same. I didn't know the author or title, but in 1990 I found copy --- it's "The Directory of Tunes" by Denys Parsons --- whilst browsing bookshelves. I excitedly bought it, thinking all my troubles were over (provided that the tune had been written before 1975, the book's publication date, of course).

How does it work? Well, consider a well known song, like "Do You Know the Way to San Jose".

Write a star for the first note: "*". Then, since the next note goes down, write a D: "*D". The third note goes up, so that's a U: "*DU". Keep going: "*DUUUD", etc. Then look up the resulting string in the book:

Perfect! Problem solved.

Except for one tiny thing. I can hear the tune in my head. But I can't tell if the notes go up, down, or stay the same (unless they change a lot). So I can't construct the string! The book has sat, unused, on my shelf for the last 20 years...

Last week I was round at some friends, and talk moved to discussing the Web, as it does. One of them was saying how much easier it was now that so many things were available to be looked up. I mentioned an anecdote told by Bertrand Meyer in 1999, about using the web to identify an opera he was listening to, and said that would only help if you have lyrics (and it would help me only if they were in English, and if I didn't fall foul of a Mondegreen).

The daughter of the house looked at me pityingly, and told me about Shazam. You point your phone in the direction of the music and "tag" it; the app records for about 10 seconds, sends it off, and the answer comes back. Perfect!

So, naturally, I downloaded the app then and there (well, after having a brisk discussion about whether it was Superman or Captain Marvel who said Shazam!), and we all spent the next few minutes playing random bits of music at it, and seeing what it could identify (despite the background noise of an excited budgie). At first I wondered if it was using a similar system to the book, but quickly realised it had to be quite different: you don't have to start at the beginning of the tune, and it not only tells you the song title, but also the artist -- and in the case of a piece of Mozart, the orchestra and conductor. So it must be matching against the actual recording. How does it work?

I found an article in the August 2006 issue of CACM which gives a brief explanation -- more technical detail can be found following the links from the wikipedia article. Essentially it looks for spectrogram peaks, takes adjacent pairs of these, does some hashing to increase the entropy, and matches the results against the music database. Many of these peaks are just noise, and so don't match. But enough do, with the added constraint that different pairs have to match at the right time intervals, to get a high quality matching system. So, a combination of a really clever algorithm and a massive database give a fantastic ability to match tunes.

But it's matching, not "recognition" as such. So it doesn't work with live music, including amateur singing (and I don't consider the noises I make to be singing as such). It's not quite the perfect system. But it's still mind-bogglingly amazingly useful.

Read More

full circle

On our recent Lake District holiday, I took a series of photos encompassing the Castlerigg stone circle. I wasn't intending to do anything with them; I was just testing out my new camera. However, after I got home, I wondered about trying to stitch them all together into a single panorama. A bit of googling later, I downloaded Hugin.

It took me a few goes to get things working -- mostly, I think, because I hadn't been thinking of taking a panorama originally, so had moved a bit between some of the shots. Still, after playing around with control points, and masking out the people, and removing the shots with the most parallax, I manage to get a panorama that includes nearly all the stones:

These are the actual photos that Hugin stitched together:

The panorama has a slightly wobbly horizon near the secondary circle, and there are some artefacts in the stitching. But given the poor quality of the original photos for this sort of thing, the result is amazingly good. And the ability to mask out certain areas (here, people) is really useful.

And all this with open source software that I downloaded of the web with the click of a button!

Read More

modelling reality

I first came across the difficulty of modelling reality when I designed a database for my books. This was back on the early 1980s when I was a physicist, and had barely even heard the word "database" before I got ViewStore for my BBC Micro. What's a book? Easy, you might think. But what about Ace Doubles? (two novels in one binding). Or a novel in three volumes? Or the same book with a variant title? Or a new edition? Or a translation? Or an abridgement? Or a book published over several editions of a magazine? Or...

The next big realisation was alphabets. In the 1990s I was involved in the international standardisation of the Z formal specification language. We needed to standardise the character set, which included lots of mathematical symbols. And our Japanese collegues wanted non-Latin alphabets. How hard could it be to allow different alphabets? So I went off to read the Unicode Standard. Oh.

I was remined of these issues when I saw a post by Charlie Stross discussing the post "Falsehoods Programmers Believe About Names" by Patrick McKenzie. People's names are difficult to model.

So is time. Different calendars. Calenders with negative and positive years, but no year zero (which is why some people ended up celebrating the 2500th anniversary of the Battle of Marathon a year too early). Changing calenders from Julian to Gregorian -- give us back our 11 days (of tax payments, that is). Changing calenders at different times in different countries. Leap years. The algorithm for leap years (and the difference between the Julian and Gregorian algorithms). The argument about whether the year 2000 should be a leap year. Time Zones. Changing time zones. Summer time (aka daylight savings time). Summer time coming and going at different days in different countries. Double summer time. Leap seconds. And so on.

One of the comments in Charlie's post references "Gay marriage: the database engineering perspective", a post with an interesting analysis of marriage database design (ignoring the name problem), and how some designs make changes harder than they need to be.

Notice the continual use of "ids" in that post. Names (even if we decide how to model them adequately) are not unique, and so are not suitable for an identifier. What properties might make a good unique identifier? I recall hearing of a new police database that used "Surname, initial, date of birth" as a unique id. The story goes that this database was installed a few days before the Kray twins were arrested...

In fact, no attribute is immutable, and so should not be used in such an identifier. The same person who told me the Kary twins story also told me of the problems a hospital had in assuming that the "sex" entry was immutable when they did their first sex-change operation.

Wait a minute. No attribute is immutable? What about "date of birth" (dob)? How can that change? Well, remember this database information is a model of reality, not reality itself. Models can have errors. The dob might have been mis-entered, or been lied about (when my maternal grandmother died, the family discovered from her birth certificate that she was 10 years older than she had let on to her children), or it might simply be unknown.

So ... if modelling designed (social) reality is this hard, why are we surprised that it is hard to model the natural world? If our own names, dates, whatever, refuse to fit some neat classification system, why should biology, geology, astronomy? Is a virus alive or not? Is something a separate species or not? Is Pluto a planet or not? When reality doesn't fit our classification, the fault lies not with reality, but with the classification system.

Read More

visualising algorithms through dance

I'm nowhere near the first to post about this, but I have to advertise this. It is a wonderful set of folk dances, from AlgoRythmics, performed by a dance troupe in association with Sapientia University, Romania, that animate six different sorting algorithms. Sheer genius.

They have produced YouTube videos showing: bubble-sort via Hungarian folk dance, insert sort via Romanian folk dance, select sort via Gypsy folk dance, shell sort via another Hungarian folk dance, merge sort via German folk dance, and, almost unbelievably, quick sort with yet another Hungarian Folk dance (using hats as index pointers!)

Read More

Total Eclipse of the Flowchart

I'm a Meatloaf/Bonnie Tyler/Jim-Steinman-as-lyricist fan. And the song "Total Eclipse of the Heart" has extra sfnal resonances for me, too. Several years ago, at a Science Fiction convention, I watched Dave Lalley's "Miscellaneous 2" tape, a compilation of 1970s and 1980s TV clips of adverts and songs. Near the end is a wonderful edit of a load of Dr Who (ToS) clips to Bonnie Tyler's TEotH -- marvellous. Nowadays there are several different YouTube videos, featuring clips of the new Doctors.

So when my other half said "Google Total Eclipse of the Heart flowchart", I did. And found the wonderful

at jeannr's blog. (Although the professional in me is obliged to point out that it isn't a flowchart, but a railroad track syntax diagram.)

Digging further through this marvellous blog, I found that the idea had been further elaborated by her followers. There's a more detailed chart

and even a video of the chart being drawn dynamically as the song is sung (by the cast of Glee, but, hey).

I love the Web.

Read More