I have found a new hobby: simulating evolutionary models. Mostly I do it in the train on my way to and from work. Yes, you can do a lot in 1 hour a day! And it’s quite captivating.

When I tell people about this, they ask me how I make simulations. Here’s how.

Yes, I program them myself. I use Octave, a free and open-source environment similar to Matlab. Octave supports a large part of the Matlab syntax and library. It also has its own extensions, such as C-like operators. I like to use ++ and +=, but today’s example will be Matlab-compatible.

The example simulates the survival of the fittest, and is only 40-odd lines long.

Here we go:

worldCapacity = 100; survivalScores = [1 1.2]; % of the 2 types initNum = worldCapacity; init2Ratio = 0.1; % The ratio of 2's initially meanNumChildren = 1; maxNumChildren = 3;

Our world has a limited capacity — more about that later. There will be two kinds of cells, type 1 and type 2. As you see, type 2 is slightly better at survival.

First, we create our population:

% init population = []; for i = 1 : initNum population(i) = (rand() < init2Ratio) + 1; end

`rand()`

returns a random number between 0 and 1. `rand() < p`

is true with probability p. `rand(i, j)`

returns an i x j matrix.

As you see, initially, only 10% of the cells are type 2.

generation = 0; ratios = []; while 1 generation = generation + 1; % show ratios(generation) = length(find(population == 2)) / length(population); plot(1 : generation, ratios); drawnow if kbhit(1) break; end

I use `kbhit()`

because Octave doesn’t always handle Ctrl-C well.

%reproduce children = []; for cell = population for i = 1 : maxNumChildren if rand() < meanNumChildren / maxNumChildren children(end + 1) = cell; end end end population = [population, children];

As you see, reproduction is asexual. With maxNumChildren large enough, the distribution becomes close to Gaussian, but that’s not today’s point.

%survive scores = survivalScores(population); prob = scores * (worldCapacity / length(population)); survive = rand(1, length(population)) < prob; population = population(survive); end

Our world has limited resources and can host only about `worldCapacity`

cells. So they have to compete for survival. There is no absolute survival in this model, i.e., if a cell is in the world alone, it won’t die. The reason a cell *can* die is that someone else is better at getting resources from the environment.

Most of the time I don’t use absolute survival at all. The reason to use competitive survival rather than (only) absolute survival is that with absolute survival, unless you calibrate the probabilities very well, you population will either die out or explode, making your simulation slow and eating up too much memory. The other reason of course is that in the real world, resources are indeed limited.

Competitive survival works like this: Each organism gets a survival score. They are then all multiplied by the same multiplier, such that their sum becomes `worldCapacity`

. The products are survival probabilities. Of course, if a “probability” is greater than 1, the organism surely survives. This may be a little unnaturalistic and also gives the population a bias towards an underpopulated world.

What next? The next step is to download the whole program and run the simulation!

July 22, 2010 at 9:33 pm

These evolutionary models you are playing with are simple models of ALife. My Thesis is in the field of Evolutionary Robotics / Multi Objective Optimization / Competitive Coevolution. We can share some ideas if you want.

July 23, 2010 at 10:55 pm

Never knew about ALife, thanks. My inspiration was “The Selfish Gene”.

Did you actually build any robots for the thesis?

January 15, 2011 at 10:51 pm

Hi Lev,

I just entered WorldPress.com for another blog I’m watching and saw that you left me a comment a long time ago.

The answer is no. I don’t build phisical robots but my simulated robots subscribe to the laws of physics.