Simulating a Tragedy of the Commons I

Tags

Code, simulations, tragedy of the commons

I introduced a recent paper, in the previous post, on the Tragedy of the Commons. The paper is published in Sustainability, an open access paper and is therefore available to everyone. One bit that I did not include in the paper, however, is the code that I used for the simulations therein. On the suggestion of a friend, I’ll publish them here in a series of posts. (Thanks Mauro)

The simulations are very simple, being based on difference equations, and were all coded in Octave. The code is not pretty, but I cleaned it up and added some comments. If you have questions, just drop me a comment on the blog. This first simulation is of a basic tragedy involving two users, and corresponds to Figures 1B-C in the paper. A resource is simulated to a steady state based on a Ricker model, for 100 steps, and then users begin to utilize it. They increase their utilization per time step at a steady rate (acceleration of utilization is constant), and the simulation is run until the resource is exhausted. The basic output of resource, one user’s benefit, and total utilization, are plotted in the figure here. The equations simulated are outlined in the excerpt from pg. 754 of the paper:
“We now introduce a term to represent consumption or utilization by our human TOC agents.

where U(t) is the total fraction of R utilized by human users at time t, and ranges from 0 to 1. The term in square brackets on the right hand side of the equation is the modified growth rate of R. The resource is stable when this term is positive or zero, that is, the rate of resource renewal exceeds or is equal to utilization. Since U is the total standardized utilization rate of all users, it may be expanded to

where there are N users and is the standardized utilization rate of the user at time t. The benefit gained from utilization is modelled as

where is the benefit to user i, f is a factor that converts resources utilized into another commodity, for example converting harvested food to energy or currency, and B is the total benefit to all users.”

The code follows, and can be executed by running it in an interactive Octave session. Please be aware that WordPress has somewhat limited options for formatting code. I’ve used Matlab formatting here, but wrapped lines might not be obvious. Single lines are all terminated with a semicolon, “;”, so just look out for those.

#RESOURCE DATA
#initial resource level
R0 = 1;
R = zeros(1,2);
R(1, 1) = 1;
R(1, 2) = 1;
#carrying capacity
K = 50;
#intrinsic growth rate r
r = 1.5;

#USER DATA
#There are 2 users
#initial acceleration of utilization
du0 = 1.05;
#individual user parameters, users 1 and 2
#initialize arrays
u1 = zeros(1,2);
u2 = zeros(1,2);
b1 = zeros(1,2);
b2 = zeros(1,2);
u1(1,1) = 1;
#initial utilization rate
u1(1,2) = .001;
#conversion factor of resource to benefit
f1 = 1;
b1(1,1) = 1;
b1(1,2) = 0;
u2(1,1) = 1;
#initial utilization rate
u2(1,2) = .001;
#conversion factor of resource to benefit
f2 = 1;
b2(1,1) = 1;
b2(1,2) = 0;
U = zeros(1,2);
U(1,1) = 1;
#total utilization
U(1,2) = b1(1,2) + b2(1,2);
du = zeros(1,2);
du(1,1) = 0;
du(1,2) = du0;

#BEGIN SIMULATION

#RUN RESOURCE TO STEADY STATE, NO UTILIZATION
for count2 = 2:100
  R(count2, 1) = count2;
  R(count2, 2) = R(count2-1, 2) * (exp(r*(1-(R(count2-1, 2)/K))));
  u1(count2, 1) = count2;
  u1(count2, 2) = u1(1,2);
  u2(count2, 1) = count2;
  u2(count2, 2) = u2(1,2);
  du(count2,1) = count2;
  du(count2,2) = du0;
endfor

#INITIATE UTILIZATION
for count1 = 101:265
  R(count1, 1) = count1;
  #modified Ricker model with utilization
  R(count1, 2) = R(count1-1, 2) * (exp(r*(1-(R(count1-1, 2)/K))) - (u1(count1-1,2)+u2(count1-1,2)));
  u1(count1, 1) = count1;
  #increase utilization
  u1(count1, 2) = u1(count1-1, 2) * du(count1-1,2);
  #benefits (equal for both users in this case)
  b1(count1, 1) = count1;
  b1(count1, 2) = f1 * u1(count1-1, 2) * R(count1-1, 2);
  b2(count1, 1) = count1;
  b2(count1, 2) = f2 * u2(count1-1, 2) * R(count1-1, 2);
  u2(count1, 1) = count1;
  u2(count1, 2) = u2(count1-1, 2) * du(count1-1,2);
  #total benefits
  U(count1, 1) = count1;
  U(count1, 2) = b1(count1, 2) + b2(count1, 2);
  du(count1,1) = count1;
  du(count1,2) = du(count1-1,2);
endfor

Python in Biodiversity Informatics

Tags

bioinformatics, Code, software

Python in Biodiversity Informatics. Nice collection of Python-based software for bio-geeks.

C++ strcpy, etc.

Tags

c++, Code, strcpy

In making some minor modifications to the main simulation program, CEG-1.4.8, I encountered a compile error using gcc 4.3.2-7. The system is 64-bit AMD based, running Fedora 10. gcc was installed from the rpm package gcc-4.3.2-7.x86_64. The simple compile command g++ -o CEG-1.4.8-SMP CEG-1.4.8-SMP.cpp yielded the following output:
CEG-1.4.8-SMP.cpp: In function ‘int main(int, char**)’:
CEG-1.4.8-SMP.cpp:75: error: ‘strcpy’ was not declared in this scope
CEG-1.4.8-SMP.cpp:114: error: ‘strcpy’ was not declared in this scope
CEG-1.4.8-SMP.cpp:281: error: ‘strcpy’ was not declared in this scope
Given that these lines refer to fairly old parts of the code, and that compilation has been successful for years, I suspected that either there is a bug in this version of gcc, this latest RH update to its gcc package, or that my code is no longer legal. It seems to be the latter, and using strcpy is not really a good idea anyway. Here’s an example of the offending code (all uses referred to file reads from command line arguments or opening files for output):
//open existing matrix file specified by user
do
{
strcpy(matrix_name,argv[2]);
std::cout << "\nName of network matrix file: " << matrix_name;
matrix_file.open(matrix_name);
}
Note the use of strcpy in the fourth line. Simply replacing the routine worked. Here’s a correction:
//open existing matrix file specified by user
do
{
const char *matrix_name = argv[2];
std::cout << "\nName of network matrix file: " << matrix_name;
matrix_file.open(matrix_name);
}

Generating beta variates

Tags

beta distribution, Code

100 Beta variates from example in text.

Our problems are to (1) draw random variates from various beta distributions, (2) embed the code in the “parallel” SMP version of CEG, and (3) pass the beta variates to CEG simulations. After exploring a number of different ways to generate the variates in C++, notably using the GNU Scientific Library, I’ve decided to simply use an Octave script. The syntax is very easy, and because the C++ simulations are embedded in a Bash script, it’s really easy to simply generate the variates and have C++ load the output files. Here’s the basic Octave script:

#! /bin/octave -qf
#Generate beta variates;
arg_list = argv ();
a = str2double(arg_list{1});
b = str2double(arg_list{2});
c = str2double(arg_list{3});
x = betarnd(a,b,c,1);
disp (x)
This script is run by the following Bash script:
#! /bin/bash
beta_file1=$4
echo $beta_file1
octave -q $1 $2 $3 >$beta_file1
The Bash script is called with command line parameters, namely the parameters for Octave’s betarnd function (a, b, and number of variates; the last parameter is 1, since output is a vector). A call would therefore look like this:
./beta.sh 2 5 100 beta_output1.dat
where beta.sh is the Bash script, the distribution being called is (2,5), and 100 random variates will be generated. The output is re-directed to output file beta_output1.dat. Don’t forget to make your script executable with chmod +x beta.sh