Optimizing a Prey-Predator model using Genetic Algorithms

Picture source: The Royal Society, Flickr

The programming section of this blog has been quite inactive lately. Hence, here’s an interesting little project I did last year.

The idea started off as a biophysics class presentation in collaboration with Abhilash Sahoo. From there, I adopted it into an optimization tool through Genetic Algorithms.

I’ve used MATLAB for all purposes. (Click here to view the code.)

The final animation:

The herding concept is inspired from Strogatz’s TED Talk: “How things in nature tend to sync up”.

Here is the basic skeleton:

  1. Green particles are prey
  2. Red particles are predators
  3. Each prey has 3 forces acting on it:
    1. a weak long-distance attractive force towards other prey
    2. a strong short-distance repulsive force away from other prey
    3. a strong repulsive force from predators
  4. The first two forces make the prey form a herd, while maintaining an equilibrium distance between themselves within the herd. The third force makes them flee from incoming predators.
  5. Each predator has 2 forces acting on it:
    1. A long-distance medium attractive force towards the general direction of a herd of prey
    2. A short-distance targeted strong attractive force towards a particular prey
  6. The predator starts its hunt using the first force. When it gets close to a prey, it switches to the second force to target the particular prey while ignoring the rest of the herd.
  7. The targeting is denoted in the video through red lines connecting the prey to the predator that appear when a prey is targeted.
  8. The green enclosure contains food for the prey.
  9. Predators have higher top speeds but are more massive. They can catch up to the prey in a straight line, but can’t turn as fast.
  10. The prey swerve automatically due to the nature of the forces, which sometimes allow them to escape.

Things get more interesting when we introduce natural selection. This can be done by creating survival parameters and the concepts of death, reproduction and imperfect copying.

Death happens when:

  1. A particle starves, that is:
    1. If a prey stays outside the green boundary for a long time
    2. If a predator cannot catch a prey for a long time
  2. A particle gets old, i.e. it runs out of its natural lifespan
  3. A prey gets caught by a predator

The survival parameters have different values for each prey/predator. I call them survival parameters as they are directly linked to the chances/duration of survival in our little virtual reality. The ones that I’ve included in this model are:

  1. The ability to accelerate (or inverse mass)
  2. Top speed
  3. Range of vision (e.g. from how far a prey can detect an approaching predator)
  4. Maximum duration of survival without food
  5. Natural lifespan
  6. The frequency of reproduction
  7. The strength of herding tendency (in case of preys)

Reproduction here is defined as a particle spawning a copy of itself. Each particle reproduces at a certain frequency. The child inherits the same values of parameters from the parent, but with small random errors. Over time these errors lead to a natural selection of parameter values in a way that maximizes the probable survival duration of the particles (since more survival duration means more number of offspring).

Leaving the model here would not be of any real use though. It’s a no-brainer that all the survival parameters will keep increasing in value over time; we don’t need a code to tell us that!

Hence we introduce interdependencies between the parameters to make them competitive. For example, we decide that a greater ability to accelerate (or less mass) would imply that the particle cannot survive without food for a long time. Thus we have created a trade-off between two desirable parameters.

This simple tweak makes this model a powerful tool of optimization. The code will ultimately maximize the survival duration and we can find out the values of the parameters that comprise the optimal solution. Basically, we will be able to quantify exactly how important each of these abstract survival parameters are for the particles in the simulation.

We can use similar techniques to optimize highly complex systems with several interdependent variables.


Blowing up a moving plane

I felt like sharing the solution of a simple yet interesting problem from the end-sem exam of our Numerical Analysis course. I’ve used R (learnt another language, yay!) and OriginPro. The solution involves solving differential equations and interpolation.


Edit: The position values in x,y,z have been given in kilometres.

Let’s break down the problem:

1. We have been given velocity equations of the plane, hence we can integrate them using the given initial position (0, 0, 0) and get the trajectory of the plane till the end of time. We can determine exactly where the plane is after 5 minutes in x, y, z coordinates. This is the time when it crosses over to the enemy country and becomes visible to the enemy radar.

2. Now we switch over to the enemy observer. He can observe the plane in polar coordinates (r, θ1, θ2). We only need to worry about r, as the missile launcher can calculate angles by itself. So we begin taking down observations for r each second. Although the value of r is readily available to the observer, we need to calculate it using the instantaneous position of the plane and position of the observer (simple distance between two points).

3. We wait till the plane enters into the sphere of 500m radius around the observer. Now the observer has to calculate where the plane will be 1 second later, and fire his missile accordingly with suitable speed. (Note that we know exactly where the plane is going, but the observer doesn’t. We cannot use our knowledge of the plane’s equation of motion.)

4. The observer needs to use the data he recorded in the past few seconds and predict a point where the plane is expected to be 1 second later.


We use Runge-Kutta algorithm to integrate the velocity equations and obtain the values of position at each second. We start tracking the distance between plane and observer after 5 mins, till the distance becomes 500m.


f <- function(t, s)
s1 <- 0.2*exp(-0.01*(t-220))/(1+exp(-0.01*(t-220)))^2
s2 <- 0.01994711*exp(-0.5*((t-200)/100)^2)
s3 <- 0.0052*exp(-t/100)
return(c(s1, s2, s3))

s0 <- c(0,0,0)
a <- 0; b <- 500
n <- b
sol <- rk4sys(f, a, b, s0, n)

time <- 301
radar <- matrix(nrow=10,ncol=2)
r <- 1000

r = sqrt((sol$y[time,1]-12)^2+(sol$y[time,2]-4.15)^2+(sol$y[time,3]^2))
r = r*1000
time = time + 1
radar[time-301,1] <- time-302
radar[time-301,2] <- r

The following file contains the calculated values of position (x, y, z) at each second starting from 0 seconds: Planepos.txt

Calculated values of r (time starting after 5 minutes):

0 s ——– 534.4396 m
1 s ——– 519.4016 m
2 s ——– 507.8468 m
3 s ——– 499.9667 m

Now we need to find a polynomial function of order 3 which fits these 4 points. It can be done using various algorithms for numerical interpolation. We then use the polynomial to find a 5th point, which will be the approximate position of the plane after 1 second (i.e. the 4th second).

Here I’ve used OriginPro to directly find a fitting polynomial after plotting the above points. However, this is just an approximate fit. The accuracy can be greatly improved with a code for interpolation.

Cubic Fit

The equation turns out to be r = 539.784 – 26.51383 t + 6.99025 t2 – 0.85882 t3.

Using the above equation, the calculated value of r at the 4th second = 490.6082 m.

Therefore, the missile needs to be fired after tracking the plane for 3 s, with a speed of 490.6082 m/s.


Probability distribution on Monopoly squares (using Python)


It is important to realize that all squares on a Monopoly board aren’t equally valuable in terms of payoff. It is common knowledge that the green and the deep blue squares are slightly less visited, whereas, red squares seem to be a lot busier. But these intuitions have no solid proof.

Here’s a program which calculates the probability distribution of a person landing on each square on a Classic Monopoly (UK) board at the end of his turns. After many faulty algorithms, I came up with what I think is a reasonably accurate version. This took approximately 10 hours of coding, pondering, and debugging. 10,000,000 data have been taken to generate probabilities.

Click to expand:

from random import *

print "This program will calculate the probabilities of each square on a Classic Monopoly (UK) board.\nThe accuracy depends on the number of calculations made.\nA large number of calculations will yield an accurate result, but will take time.\n"

matrix = [[0 for i in range(40)] for j in range(40)]
square_names=["Go", "Old Kent Road", "Community Chest", "Whitechapel Road", "Income Tax", "King's Cross station", "The Angel Islington", "Chance", "Euston Road", "Pentonville Road", "Jail", "Pall Mall", "Electric Company", "Whitehall", "Northumberland Avenue", "Marylebone station", "Bow Street", "Community Chest", "Marlborough Street", "Vine Street", "Free Parking", "Strand", "Chance", "Fleet Street", "Trafalgar Square", "Fenchurch Street station", "Leicester Square", "Coventry Street", "Water Works", "Piccadilly", "Go to Jail", "Regent Street", "Oxford Street", "Community Chest", "Bond Street", "Liverpool Street station", "Chance", "Park Lane", "Super Tax", "Mayfair"]
turns=input("Enter number of calculations: ")

def setdice(n):
    print "\nSetting dice probabilities..."
    f1 = open("dice_probabilities.dat", 'w')
    for i in range(n):
        if a==b:
            if c==d:
                if e!=f:
                    cnt = a+b+c+d+e+f
    for i in range(36):
    print "Dice probabilities have been set using Monopoly rules.\n"+str(turns)+" steps have been used.\nData saved to dice_probabilities.dat\n"


for i in range(40):
    for j in range(40):
        if k<0:
        if k<36:

def chance(x):
    if n==1:
        return 39
    elif n==2:
        return 0
    elif n==3:
        return 30
    elif n==4:
        return x-3
    elif n==5:
        return 24
    elif n==6:
        return 15
    elif n==7:
        return 11
        return x

def community_chest(x):
    if n==1:
        return 0
    elif n==2:
        return 30
    elif n==3:
        return 1
        return x

def square_probs(num):
    print "Calculating square probabilities..."
    while x<num: n=randint(1,num) i=0 for i in range(40): b=a+i if b>39:
            if n<=0: break a=b if a in [7,22,36]: a=chance(a) if a in [2,17,33]: a=community_chest(a) squares[a] = squares[a]+1 x=x+1 if a==30: a=10 squares[30]=0 squares[10]=squares[10]+2 x=x+1 f2 = open("square_probabilities.dat",'w') for i in range(40): probability[i]=squares[i]/float(num) print>>f2,(str(i+1)+"\t"+str(probability[i])+"\t"+square_names[i])
    print "Square probabilities have been calculated using " +str(num)+" dice rolls.\nData saved to square_probabilities.dat"

This program will calculate the probabilities of each square on a Classic Monopoly (UK) board.
The accuracy depends on the number of calculations made.
A large number of calculations will yield an accurate result, but will take time.

Enter number of calculations: 10000000

Setting dice probabilities...
Dice probabilities have been set using Monopoly rules.
10000000 steps have been used.
Data saved to dice_probabilities.dat

Calculating square probabilities...
Square probabilities have been calculated using 10000000 dice rolls.
Data saved to square_probabilities.dat

Approximations made:

  1. The program considers that a player rolls only once per turn. In reality, a player gets to roll again in case of doubles. To minimize error, I’ve come up with something which I’ll call a pseudo-multi-roll system. It will be explained shortly. This method is a leftover from one of my previous algorithms. I kept it because I believe it is faster than actual multiple rolls.
  2. A player does not buy his way out of jail, i.e. he stays in jail for three turns whenever he gets sent to it. Thus the probability of the jail square calculated is more than its actual value.
  3. Chance and Community Chest cards are completely reshuffled after each draw. This increases the probability of squares which have associated cards (e.g. Go, Mayfair, Trafalgar Square, Pall Mall, Old Kent Road, Jail etc.). This assumption is an error only if the drawn cards are (lazily) kept at the bottom of the decks in actual games.
  4. The probability is calculated by considering the position of the player at the end of each turn. Hence, the square “Go to Jail” has 0 probability, since no one occupies it after the turn is over. Also, the probability of the Chance and Community Chest squares decreases a lot because they have cards sending players to other squares.

The Pseudo-Multi-Roll system:

The system takes into account that a player rolls again if he gets a double (and that he can roll a maximum of two doubles per turn), and generates a probability distribution after considering 10,000,000 rolls. This probability distribution gives us the chance of getting a particular sum at the end of a turn. For example, we can know the probability of the player getting a sum of 17 (say the dice gave [6,6] on first roll and [2,3] on second). We will not know that he had stopped at 12 in the middle, before advancing to 17, however this should not affect the end result considerably. The system lowers the probability of “Jail” since it doesn’t consider triple doubles.

The data from “dice_probabilities.dat” has been plotted. You can observe a slight increase of probabilities for the numbers between 7 and 12 due to the pseudo-multi-roll.

monopoly dice

The Probability Matrix:

This is a 40×40 table, each row of which will store a copy of the probability chart generated above. It is stored in a way that the element at the j-th column of the i-th row will give the probability of travelling from square number ‘i’ to square number ‘j’ in one turn. This matrix will be referred after each turn to determine where the player will go next (by weighted random generation method).

Fine Tuning:

“Chance”, “Community Chest” and “Go to Jail” have been hard-coded separately to generate their desired results.

Probabilities of each square:

This has been done by playing a game of 10,000,000 turns. The number of times the player stays on each square at the end of his turns has been recorded, and the probabilities have been calculated accordingly.

The data from “square_probabilities.dat” has been plotted.


Following observations are crucial:

  1. The red and orange squares are indeed the best squares on the board in terms of probability of visits. This is due to their location after jail (which draws in players from other parts of the board). The pink squares don’t get that big an advantage from jail, since they’re located too close (the dice probabilities are highest for values between 4 and 10).
  2. Jail has a drastically high probability because there are many ways of landing in jail. Also, one landing consumes three turns.
  3. The Chance and Community Chest squares show less probability since they can send the player to different squares within the same turn.
  4. The Go to Jail has 0 probability since it sends the player to jail.
  5. The squares with Chance/Community Chest cards have slightly higher probabilities (e.g. Go, Mayfair, Trafalgar Square, Pall Mall, Old Kent Road etc.).
  6. The squares before “Jail” and after “Go to Jail” have slightly lower probabilities because they’re skipped whenever the player lands on “Go to Jail”.
  7. “Marylebone station” is the best station and “Water Works” is the better utility.

The Glowing Bulb Problem

I came across this interesting problem yesterday.

There are few light bulbs and all of them are switched off. On the first turn, you press the switch of all bulbs, i.e. you turn all of them on. On the second turn, you switch all the second bulbs (2, 4, 6…). On the third turn, you switch all third bulbs, and so on.

An example with five bulbs and five steps (where 1 denotes a glowing bulb, and 0 denotes a non-glowing bulb):

00000 (initial)
11111 (step 1)
10101 (step 2)
10001 (step 3)
10011 (step 4)
10010 (step 5)

Hence, in case of five bulbs, the number of glowing bulbs at the end is 2. You have to find out how many bulbs are left glowing after 100 turns when you have 100 bulbs. The solution must be done mentally, without any rigorous calculation.

It is recommended that you try it yourself before reading the solution.

Mental Solution:

Every bulb is off at the beginning. It is apparent that it will remain off if it is switched even number of times. Likewise, it’ll be on if switched odd number of times. Let’s pick a random bulb, say bulb 8. Bulb 8 will be switched on the 1st step, 2nd step, 4th step and the 8th step, i.e. it’ll be switched on a step if the step number is one of its factors.

Merging these two conditions, we have: a bulb will remain on if it has an odd number of factors.

Now we have the task of finding out which digits have odd number of factors. We can always divide a number into a set of two factors. For example, 8 can be divided into (1,8) and (2,4). These factors multiply with each other to give the original number. So we’ll always have even number of factors, unless we get two equal factors, i.e.  (x,x). Only perfect squares can be divided into two equal factors.

8 factors                                  16 factors

So the problem boils down to: “How many perfect squares are there under (and including) 100?” The answer is 10. 10 bulbs will remain glowing after the 100th step. (The logic works for any number of bulbs, if the number of steps is equal to the number of bulbs.)

Programmed Solution (using Python):

Such problems can be easily solved through computer programs. It leads to some interesting results.

n=input("Enter number of bulbs: ")
m=input("Enter number of steps: ")
if (m>n):
f1 = open("lights.dat", 'w')
for x in range(m):
    for i in range(x,n,x+1):
    for y in range(n):
    print "\n" * 100
    print "Executing... " +str(float((x+1))*(100.0/m))+ "% completed."
print "The file has been created."
print "Number of glowing bulbs = "+ str(cnt)


Enter number of bulbs: 100
Enter number of steps: 100

Executing… 100.0% completed.
The file has been created.
Number of glowing bulbs = 10

The data has been written on a file “lights.dat” which when plotted gives:

We can see that 10 bulbs are left glowing after the 100th step. We can also guess at the logic by analyzing the graph. The first region is quite chaotic and we ignore it. Then, there’s a region (before 49) where the number of glowing bulbs remain constant. This is the region where one bulb is turned on and one is turned off in each step. At 49, two bulbs get turned on, resulting in a spike (our first clue, what’s special about 49?). After 50 comes the steady decline, where one bulb is turned off in each step. But as you see, in the each of the steps 64, 81, and 100 one bulb is turned on instead. We can immediately guess that bulbs numbered as perfect squares get turned on at their own step number.

When I opened the data file, I found something remarkable. By the way, the columns here are:

  1. the step number
  2. number of glowing bulbs
  3. the status of bulbs (a 100 digit number containing 0s and 1s denoting the state of each bulb in order)

Lights data
The pattern of 0s and 1s is obvious. The vertical stripes (green) are obviously the perfect squares. The lower red line is formed because bulbs (except perfect squares) get turned off on their own step number. The upper red line is formed by the halves of the even numbers which turn on the respective numbers’ bulbs.

Note that every number apart from perfect squares get turned off at the lower red line. No change takes place in any bulb below that, and hence, only perfect squares stay glowing.


Simulation of Random Walk Hypothesis (using Python)

The Random Walk Hypothesis (as stated in Feynman’s Lectures) has always troubled me. Maybe the problem lies in my basic understanding of probability. Whatever it is, I can’t seem to get a grip on it.

There is a man who tosses a coin. If heads, he moves one step forward. If tails, he moves one step backward. Where will he end up after N number of tosses? Let us call his final distance as D. We are considering a sufficiently large N, i.e. N → ∞. We might think that D → 0, since eventually the number of heads and number of tails will balance out (both of their probabilities being equal). That is not the case.

For example,

when N=2, we might easily to get 2 heads and 0 tails, so D=2.
when N=50, we might get 28 heads and 22 tails, so D=6.
when N=1000, 520 heads and 480 tails, D=40.
and so on…

The above numbers are just speculations using common sense. Note that D does not decrease, but is rather likely to increase at a slow rate, with increasing N. Due to the difference in rates of increasing, the quantity D/N is likely to decrease with increasing N.

Hence, D/N → 0 when N → ∞.

The Random Walk Hypothesis states that the expectation value of D () should be √N. You can see that it agrees with the condition: D/N → 0 for N → ∞. The mathematical proof is simple too. What bothers me is that this expectation value is seldom reflected in experiments.

In the following simulation, I’ve assumed that Python’s randomness generating skills are good enough for the task. Any number of experiments can be performed. Each experiment will consist of an user-defined number of steps. The program will calculate the distance in each experiment, and display the number of experiments which matched the theory, and the absolute average distance. Ultimately it will state whether the average experiment agrees with the theory or not. I’ve allowed a standard error of 10%.

import random
noe=input("Enter number of experiments: ")
steps=input("Enter number of steps per experiment: ")
print "Performing experiments..."
for a in range(0,noe):
    for i in range(0,steps):
        if j==1:
    if apos<=(1.1*(steps**(0.5))) and apos>=(0.9*(steps**(0.5))):
print "Expected value of distance: %d"%(steps**(0.5))
print "Number of experiments agreeing with the theory = %d"%cnt
print "Absolute average distance: "+ str(avg)
if cnt>=(noe-cnt):
    print "The average experiment satisfies the theory."
    print "The average experiment disagrees with the theory."

Sample Output:

Enter number of experiments: 50
Enter number of steps per experiment: 1000000
Performing experiments…
Expected value of distance: 1000
Number of experiments agreeing with the theory = 4
Absolute average distance: 801.32
The average experiment disagrees with the theory.

As you can see, the average experiment with N = 1 million steps each failed to agree with the theory. But the average distance (801.32) looks promisingly close to 1000, so higher values of N might work. It’ll take a long time to execute though.


Riju’s Vault System (using Python)

I am a lazy person. I don’t code much. When I get an idea, I spend time pondering whether it’ll be possible to implement it in a code within the limits of my programming knowledge. If yes, then it’s as good as done. If no, then there’s nothing to be done since my ideas are never worth learning new things.

But today’s an exception. I wrote a program! Very basic, though. No groundbreaking concept, just a product of loops and objects and boredom.

It can be done in a simpler way probably.

The program creates a system of vaults. To create a vault, a person needs to enter his name and set a password. A vault can be accessed using its serial number and password. Once access is granted to a vault, one can check its status, and make deposits and withdrawals. It is necessary to sign out in order to close a vault.

Part of the code is mere decoration (time delays and stuff). I’ve checked almost all possible scenarios and found it to be error free, so do notify me in case you find some bugs.

import time
class Vault:

    def create(self,name,passw):
        print("Creating vault...")
        print("Vault number %d created in the name of %s." %(self.sno,

    def check(self):
        print("Fetching balance info...")
        print("Your current balance is Rs.%d"%self.balance)
        print " "

    def authenticate(self,passw):
        if self.password == passw:
            print("Access granted!")
            return 1
            print("Access denied! Invalid password.")
            print " "
            return 0

    def deposit(self,a):
        print("Rs.%d has been deposited to your vault. Your current balance is Rs.%d."%(a,self.balance))
        print(" ")

    def withdraw(self,a):
        if self.balance>=a:
            print("Withdrawal successful! New balance is Rs.%d." %self.balance)
            print(" ")
            print("Insufficient balance.")
            print(" ")

print "Hello! Welcome to Riju's vault system."
while 1:
    print "1. Create a vault"
    print "2. Sign in"
    print "3. Exit"
    choice = input("Enter choice: ")
    if choice == 1:
        v = Vault()
        tname = raw_input("Enter your name: ")
        tpass = raw_input("Set your password: ")
        tamount = input("Enter starting deposit: ")
    elif choice == 2:
        tsno = input("Enter vault number: ")
        if len(Vault.vaultlist)<tsno:
            print("Vault does not exist.")
            print " "
        tpass = raw_input("Enter password: ")
        if Vault.vaultlist[(tsno-1)].authenticate(tpass) == 1:
            print("Welcome %s." %Vault.vaultlist[(tsno-1)].name)
            while 1:
                print("1. Check balance")
                print("2. Deposit")
                print("3. Withdraw")
                print("4. Sign out")
                choice = input("Enter choice: ")
                if choice == 1:
                elif choice == 2:
                    tamount=input("Enter amount to deposit: ")
                elif choice == 3:
                    tamount=input("Enter amount to withdraw: ")
                elif choice == 4:
                    print "Signing out..."
                    print "You have been successfully signed out."
                    print " "
                    print "Invalid choice"
                    print " "
    elif choice == 3:
        print "Exiting..."
        print("Thank you for choosing Riju's vault system. You suck!")
        print "Invalid choice"

Sample output:

Hello! Welcome to Riju’s vault system.
1. Create a vault
2. Sign in
3. Exit
Enter choice: 1
Enter your name: Debapriyo
Set your password: iloveghibli
Creating vault…
Vault number 1 created in the name of Debapriyo.
Enter starting deposit: 9000
Rs.9000 has been deposited to your vault. Your current balance is Rs.9000.

1. Create a vault
2. Sign in
3. Exit
Enter choice: 2
Enter vault number: 1
Enter password: iloveghibli
Access granted!
Welcome Debapriyo.
1. Check balance
2. Deposit
3. Withdraw
4. Sign out
Enter choice: 1
Fetching balance info…
Your current balance is Rs.9000

1. Check balance
2. Deposit
3. Withdraw
4. Sign out
Enter choice: 2
Enter amount to deposit: 2000
Rs.2000 has been deposited to your vault. Your current balance is Rs.11000.

1. Check balance
2. Deposit
3. Withdraw
4. Sign out
Enter choice: 3
Enter amount to withdraw: 5000
Withdrawal successful! New balance is Rs.6000.

1. Check balance
2. Deposit
3. Withdraw
4. Sign out
Enter choice: 4
Signing out…
You have been successfully signed out.

1. Create a vault
2. Sign in
3. Exit
Enter choice: 3
Thank you for choosing Riju’s vault system. You suck!