In [529]:
%matplotlib inline

Imports


In [530]:
import random
from collections import Counter
import matplotlib.pyplot as plt
import numpy as np

Define classes


In [531]:
class resource:
    """
    """
    def __init__(self, name, amount = 10):
        self.name = name
        self.quantity = amount
    def add(self, amount):
        self.quantity += amount
    def consume(self, amount):
        self.quantity -= amount
    def get_value(self):
        return self.quantity
    
class agent:
    """
    """
    def __init__(self, name, require, output):
        if len(require) > len(set(require)):
            raise ValueError("duplicates in require list")
        if len(output) > len(set(output)):
            raise ValueError("duplicates in output list")
        self.name = name
        self.r = require
        self.o = output
    def requires(self):
        return self.r
    def outputs(self):
        return self.o

Convience functions


In [532]:
def constrained_sum_sample_pos(n, total):
    """Return a randomly chosen list of n positive integers summing to total.
    Each such list is equally likely to occur."""

    dividers = sorted(random.sample(xrange(1, total), n - 1))
    return [a - b for a, b in zip(dividers + [total], [0] + dividers)]    
    
def print_community(pop):
    for k,v in Counter(pop).items():
        print k.name, v
        
def replinish(initial_resources):
    for r in initial_resources:
        r.add(r.get_value())

Set parameters


In [533]:
time = 5000

#total_resources = 1000000
#x = constrained_sum_sample_pos(num_resources, total_resources)

#for r in num_resources:
#    resource_pool.append((name,amount))

#set resources, global
a = resource("a", 125000)
b = resource("b", 125000)
c = resource("c", 200000)
d = resource("d", 200000)
e = resource("e", 100000)
f = resource("f", 250000)
g = resource("g", 100000)
h = resource("h", 100000)

resource_pool = [a,b,c,d,e,f,g,h]
initial_pool = resource_pool[:]

#set species
s1 = agent("s1", [(a,15), (c,50), (h,5), (g,5)], [(b,10)])
s2 = agent("s2", [(b,25), (d,15), (f,5)], [(c,10)])
s3 = agent("s3", [(c,50), (g,5), (e,5)], [(d,5)])
s4 = agent("s4", [(d,15), (f,5)], [(a,20)])
s5 = agent("s5", [(a,15), (e,5)], [(b,10)])
s6 = agent("s6", [(b,25), (h,5)], [(c,20)])
s7 = agent("s7", [(c,5), (a,15), (f,5)], [(d,5)])
s8 = agent("s8", [(d,15), (h,5), (f,5)], [(a,20)])

species = [s1, s2, s3, s4, s5, s6, s7, s8]
population = []

#initialize resource counts, this is used for tracking resources over time
stepdict = {}

for i in resource_pool:
    stepdict[i.name] = [i.get_value()]

#generate abundance randomly
x = constrained_sum_sample_pos(len(species), 400)

i = 0
for s in species:
    population.extend([s for y in xrange(x[i])])
    i=i+1


#population.extend([s1 for y in xrange(x[0])])
#population.extend([s2 for y in xrange(x[1])])
#population.extend([s3 for y in xrange(x[2])])
#population.extend([s4 for y in xrange(x[3])])

random.shuffle(population)

popdict = {}
for k, v in Counter(population).items():
        popdict[k.name] = [v]

Run the model


In [534]:
print "initial population:"
print_community(population)
#run
for step in xrange(time):
    if step % 100 == 0:
        replinish(initial_pool)
        
    new_pop = []
    
    for i in population:
        required_resources = i.requires() #returns a list of tuples
        output_resources = i.outputs()
        consume = []
        double = []
        for resource, amount in required_resources:
            if resource.get_value() >= amount:
                #resource.consume(amount)
                consume.append(resource)
                if resource.get_value() >= (2 * amount):
                    double.append(resource)
        if len(double) == len(required_resources) and step % 50 == 0:
            new_pop.append(i)
        if len(consume) == len(required_resources):
            for resource, amount in required_resources:
                resource.consume(amount)
            
            for output, amnt in output_resources:
                output.add(amnt)
        else:
            population.remove(i)
    
    
    #track resources for plotting
    for i in resource_pool:
        stepdict[i.name].append(i.get_value())
    
    for k, v in Counter(population).items():
        popdict[k.name].append(v)
    
    if len(population) == 0:
        break
    population.extend(x for x in new_pop)

print "steps:", step+1, "out of", time

print "final population size:", len(population)
print_community(population)

print "final resources:"
for i in resource_pool:
    print i.name, i.get_value()


initial population:
s3 17
s8 24
s6 36
s4 136
s1 13
s7 90
s2 33
s5 51
steps: 146 out of 5000
final population size: 0
final resources:
a 5
b 15
c 35
d 5
e 124330
f 244415
g 287870
h 162550

Plots


In [535]:
labs = []
for k, v in stepdict.items():
    plt.plot(v)
    plt.title("resources")
    labs.append(k)
plt.legend(labs, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)


Out[535]:
<matplotlib.legend.Legend at 0x125b309d0>

In [536]:
labs = []
for k, v in popdict.items():
    plt.plot(v)
    plt.title("population")
    labs.append(k)
#plt.legend(labs, loc=2)
plt.legend(labs, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)


Out[536]:
<matplotlib.legend.Legend at 0x125b29690>

In [536]:

Copyright (C) 2014 Kathryn Iverson

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.


In [536]: