I've modeled the population changes for a group of agents that inhabit a 2D lattice (with wrapped boundaries) made up of M x M
number of grids. At each time step:
- A parent reproduces offsprings that disperse independently -- the number depends on the grid location.
- The offsprings live if they land on a "patch" grid but are removed otherwise.
- Neighboring offsprings compete based on local density.
- Survived offsprings grows up, becoming parents the next time step.
As an inexperienced programmer, my current Python script is pretty inefficient, taking as long as 19 minutes(!!) for the given parameters. I suspect I'm relying too much on NumPy functions but cannot see clear ways to optimize it.
import numpy as np
import random
import math as mp
import matplotlib.pyplot as plt
import time
# Lattice parameters
M = 60
patches = np.random.choice(2, M * M).reshape(M, -1)
# Biological parameters
mu = 14.5
alpha = .1
sigma = 5.
def reproduction(a_):
da = np.floor(a_) # discretize parent locations
z = np.random.normal(0, 0.1, M**2).reshape(M,-1) # generate noise in reproduction
muloc = np.ones(len(a_)) # mean reproduction rate per location
for i in xrange(len(a_)):
muloc[i] = mu * da[i,0]**0.01 + z[da[i,0], da[i,1]] # grid-dependency
if muloc[i]<0: # mean reproduction rate cannot be negative
muloc[i]=0
return muloc
def dispersal(self):
return np.random.multivariate_normal([self[0], self[1]], [[sigma**2*1, 0], [0, 1*sigma**2]])%M # Gaussian dispersal
def landed(offspring_list):
Inlist = [] # landed offspring list
og = np.zeros((M, M)) # lattice defined by the number of landed offsprings per grid
dl = np.floor(offspring_list) # discretize offspring locations
for i in range(len(offspring_list)):
if patches[dl[i,0],dl[i,1]] == 1:
Inlist.append(offspring_list[i])
return Inlist
# competition kernel
ker = np.zeros((M, M))
for i in range(M):
for j in range(M):
r2 = min((i-0), M-(i-0))**2 + min((j-0), M-(j-0))**2
rbf = mp.sqrt(1/(2*sigma**2))
ker[i,j] = np.exp(-(rbf*r2)**2)
ker = ker/np.sum(ker)
def competition(offspring_list):
dl = np.floor(offspring_list) # discretize offspring location
og = np.zeros((M,M)) # MxM lattice defined by the number of seeds at each grid
for i in range(len(offspring_list)):
og[dl[i,0],dl[i,1]] += 1
v1 = np.fft.fft2(ker) # fast fourier transform of kernel
v2 = np.fft.fft2(og) # fast fourier transform of offspring densities across lattice
v0 = np.fft.ifft2(v1*v2) # convolution via inverse fast fourier transform
dd = np.around(np.abs(v0), decimals=3) # rounding resulting values
alive = [] # indicators for competition outcome: 1 for alive, 0 for dead
density = np.ones((len(dl),1))
for i in range(len(dl)):
density[i] = dd[dl[i,0],dl[i,1]] #FFT-approximated density at each grid
mortality = 1 - 1 / (1 + density[i]/1.) # competition-based mortality rate
if np.random.random() < mortality:
alive.append(0)
else:
alive.append(1)
return alive, og, dd
# Simulation
sim_num = 5
terminal_time = 50
start_time = time.time()
for s in range(sim_num):
print 'simulation number ', s+1, ' of ', sim_num
# initial conditions
n = 200
parent = np.random.random(size=(n, 2)) * M
pop_density = [n/float(M**2)]
for t in range(terminal_time):
print 't = ', t+1
muloc = reproduction(parent)
offspring = []
for parent_id in range(len(parent)):
for m in range(np.random.poisson(muloc[parent_id])): # offsprings reproduced
offspring.append(dispersal(parent[parent_id])) # offspring dispersed
offspring = landed(offspring) # offsprings landed
if len(offspring) == 0: # possible extinction
parent = []
else:
offspring = np.vstack(offspring)
n = len(offspring)
alive = competition(offspring)[0] # offspring in competition
indexes = [i for i, x in enumerate(alive) if x==1]
if len(indexes) == 0: # possible extinction
parent = []
else:
parent = np.vstack([offspring[l] for l in indexes])
n = len(parent)
pop_density.append(n/float(M**2))
plt.plot(range(terminal_time+1), pop_density, 'k--', alpha=0.4)
print time.time() - start_time, 'seconds'
plt.xlabel('Time')
plt.ylabel('Density of Individuals')
plt.show()