Return to the article
#!/usr/local/bin/perl
# Listing_2.pl
=head1 NOTES
This script gives multiple runs (100) of the simulator, for statistical analysis. For a toy implementation
designed to run a single simulation to give the user a feel for what is going on, Listing_1.pl is the most appropriate.
As with all the scripts, Perl must be installed.
The usage is:
perl -w Listing_2.pl
Defaults will give quite a high cultural isolation.
Output will just be tab-delimited text eg:
pop size gene dispersal meme dispersal
run 0 1411 3 0
run 1 1673 2 2
run 2 1405 2 4
Redirection to a spreadsheet is necessary for analysis.
To see a rather different effect, try:
perl -w Listing_2.pl -contag 0.49 -nn_contag 0.6 -cultsel Y -natsel Y
Those parameters should cause very little meme isolation
If you want each simulation to run for longer (ie. still 100 runs but of 200 generations/run) try:
perl -w Listing_2.pl -time 200
Other things can be changed by rewriting the script, like:
$runs : the number of simulation runs in the batch (default 100)
$popsize : the initial population size (default 200)
$rep_rate : the reproduction rate (default 1%)
$capacity : the max. population capacity per cell (default quasi-infinite)
This script appears to run reasonably quickly, eg: a single run takes about 10 secs
The entire set of 100 runs will take several minutes. (Pentium III, 1GHz, running RedHat)
But toggling either -cultsel or -natsel to Y slows the simulation by a factor of 2 or so,
as natural selection increases the numbers in the population.
As with all Perl scripts, they are offered as is, without any warranties concerning
safety, functionality etc
The author welcomes any comments on any aspects of the software.
to: dgatherer2002@yahoo.co.uk
This is an occasionally on-going project and newer versions of code are frequently
available, so if you want the really rough cuts, please ask.
=cut
#------PRELIMINARIES-------------------------------------------------#
use strict;
use POSIX;
srand(); # initialise randomiser
my $start_popsize = my $popsize = 200; # initial pop
my $time = 100; # over gens
my $line = my $cols = 10; # in a 10 by 10 landscape
my $rep_rate = 0.01; # pop grows per generation
my $contag_rate = 0.14; # population contage to neighbour per gen
my $nn_contag_rate = 0.73; # contag to other culture
my $mig_rate = 0.5; # pop move every gen
my $capacity = 100000000; # max pop per cell
my $pause = 1; # when to stop and draw the map
my @pop;
my $cult_sel = "N"; # no cultural selection
my $nat_sel = "N"; # no natural selection (of trait)
for(my $x=0; $x<=$line-1; $x++)
{
for(my $y=0; $y<=$cols-1; $y++)
{
$pop[$x][$y] = 0;
}
}
foreach(@ARGV)
{
if(/-help/)
{
die "\nUsage: perl -w Listing_2.pl -num -time -size -cultsel -natsel -contag -migrate -nn_contag\n";
}
}
my %args = @ARGV;
if(exists($args{-num})) {$popsize = $args{-num}; print "\npopsize\t$popsize";}
if(exists($args{-time})) {$time = $args{-time}; print "\ntime\t$time";}
if(exists($args{-cultsel})) {$cult_sel = $args{-cultsel}; print "\ncult. sel.\t$cult_sel";}
if(exists($args{-natsel})) {$nat_sel = $args{-natsel}; print "\nnat. sel.\t$nat_sel";}
if(exists($args{-contag})) {$contag_rate = $args{-contag}; print "\ncontag.\t$contag_rate";}
if(exists($args{-nn_contag})) {$nn_contag_rate = $args{-nn_contag}; print "\nnn_contag.\t$nn_contag_rate";}
if(exists($args{-migrate})) {$mig_rate = $args{-migrate}; print "\nmigrate\t$mig_rate";}
foreach($popsize, $time)
{
if(/[\D+]/)
{
die "\npopulation size, time must all be positive integers\n";
}
}
foreach($cult_sel, $nat_sel)
{
unless(/[YN]/)
{
die "\nCultSel and NatSel must be just Y or N\n";
}
}
foreach($contag_rate, $nn_contag_rate, $mig_rate)
{
unless($_ =~ /[\d\.-]/ && ($_ <= 1) && ($_ >= 0))
{
die "\ncontagion rates and migration rate must be between 0 and 1\n";
}
}
print "\nrep rate $rep_rate";
print "\nmig rate $mig_rate";
print "\nintracell contag rate $contag_rate";
print "\nintercell contag rate $nn_contag_rate";
print "\ncultural seln $cult_sel";
print "\nnatural seln $nat_sel";
#------MAIN CODE SEQUENCE---------------------------------------------#
print "\n\tpop size\tgene dispersal\tmeme dispersal";
for(my $runs=0; $runs<=100; $runs++)
{
$popsize = $start_popsize;
my $pop = &evenly_initialise(); # or spread evenly
print "\nrun $runs";
my $g_disp = my $t_disp = 0;
for(my $gen=1;$gen<=$time;$gen++) # for a certain number of cycles
{
$pop = &reproduce($pop); # replicate and pass traits by contagion
$pop = &migrate($pop); # move strings over grid
if($gen%$pause ==0)
{
($g_disp, $t_disp) = &pop_state($pop);
} # look at it
}
print "\t$popsize";
print "\t$g_disp\t$t_disp";
}
#-----SUBROUTINES-----------------------------------------------------#
#----------------------------------------------------------------------
sub pop_state() # prints current pop
{
my $input = $_[0];
my @current_pop = @{ $input ->{POPULATION} };
my @gene_state; # arrays for calculating the
my @trait_state; # horizontal and vertical dispersal
for(my $x=0; $x<=$line-1; $x++)
{
for(my $y=0; $y<=$cols-1; $y++)
{
my $freq_A = my $freq_1 = 0;
my $cell_count = 0; # how many per cell
foreach(@current_pop)
{
my $xcoord = $_ -> { LOCATION }[0];
my $ycoord = $_ -> { LOCATION }[1];
if($xcoord == $x && $ycoord == $y)
{
$cell_count++;
if($_ -> { TRAIT } eq "A")
{
$freq_A++;
}
if($_ -> { SERIAL_NUM } == 1)
{
$freq_1++;
}
}
}
$pop[$x][$y] = $cell_count;
if($freq_A == 0 && $freq_1 == 0 && $pop[$x][$y] == 0)
{
$trait_state[$x][$y] = 0;
$gene_state[$x][$y] = 0;
}
else
{
$freq_A /= $pop[$x][$y];
$freq_1 /= $pop[$x][$y];
$trait_state[$x][$y] = $freq_A;
$gene_state[$x][$y] = $freq_1;
}
}
}
my($g, $t) = &disp_calc(\@gene_state, \@trait_state);
return ($g, $t);
}
#--------------------------------------------------------------------
sub reproduce() # rather similar to recombine
{
my $input = $_[0];
my $contag_events = my $repro_events = 0; # how many per cycles
for(my $b=0;$b<=$popsize-1;$b++) # run through pop
{
my $rand = rand();
if($nat_sel =~ /Y/i && @{ $input ->{POPULATION} }[$b] ->{TRAIT} eq "A")
{
$rand /= 2; # double the random number
}
if($rand < $rep_rate) # also duplicate every 10th one
{
my $xcoord = @{ $input ->{POPULATION} }[$b] -> { LOCATION }[0];
my $ycoord = @{ $input ->{POPULATION} }[$b] -> { LOCATION }[1];
if($pop[$xcoord][$ycoord] < $capacity)
{
my @new_one_seq;
my $new_seq =
{
LOCATION => [ $xcoord, $ycoord ],
SERIAL_NUM => @{ $input ->{POPULATION} }[$b] ->{SERIAL_NUM},
TRAIT => @{ $input ->{POPULATION} }[$b] ->{TRAIT},
};
@{ $input ->{POPULATION} }[$popsize++] = $new_seq;
$pop[$xcoord][$ycoord]++;
$repro_events++;
}
}
$rand = rand();
if($cult_sel =~ /Y/i && @{ $input ->{POPULATION} }[$b] ->{TRAIT} eq "A")
{
$rand *= 2; # halve the random number
}
if($rand < $contag_rate) # and contage every 10th one
{
my $chosen_one = floor($rand*$popsize);
while(abs(@{ $input ->{POPULATION} }[$b] -> { LOCATION }[0]
- @{ $input ->{POPULATION} }[$chosen_one] ->{ LOCATION }[0]) > 1
|| abs(@{ $input ->{POPULATION} }[$b] -> { LOCATION }[1]
- @{ $input ->{POPULATION} }[$chosen_one] -> { LOCATION }[1]) > 1)
{ # replace with random other string
my $chos_rand = rand();
$chosen_one = floor($chos_rand*$popsize);
}
if($rand < $nn_contag_rate)
{
if(abs(@{ $input ->{POPULATION} }[$b] -> { LOCATION }[0]
- @{ $input ->{POPULATION} }[$chosen_one] ->{ LOCATION }[0]) <= 1
&& abs(@{ $input ->{POPULATION} }[$b] -> { LOCATION }[1]
- @{ $input ->{POPULATION} }[$chosen_one] -> { LOCATION }[1]) <= 1)
{
@{ $input ->{POPULATION} }[$b] ->{TRAIT} = @{ $input ->{POPULATION} }[$chosen_one] ->{TRAIT};
my $xcoord = @{ $input ->{POPULATION} }[$b] -> { LOCATION }[0]; # direct access, better ??
my $ycoord = @{ $input ->{POPULATION} }[$b] -> { LOCATION }[1]; # yes, same method for string access
my $trait = @{ $input ->{POPULATION} }[$b] -> { TRAIT };
$contag_events++;
next; # only one contag per indiv per gen
}
}
else
{
while(abs(@{ $input ->{POPULATION} }[$b] -> { LOCATION }[0]
- @{ $input ->{POPULATION} }[$chosen_one] ->{ LOCATION }[0]) > 0
|| abs(@{ $input ->{POPULATION} }[$b] -> { LOCATION }[1]
- @{ $input ->{POPULATION} }[$chosen_one] -> { LOCATION }[1]) > 0)
{ # replace with random other string
my $chos_rand = rand();
$chosen_one = floor($chos_rand*$popsize);
}
@{ $input ->{POPULATION} }[$b] ->{TRAIT} = @{ $input ->{POPULATION} }[$chosen_one] ->{TRAIT};
my $xcoord = @{ $input ->{POPULATION} }[$b] -> { LOCATION }[0]; # direct access, better ??
my $ycoord = @{ $input ->{POPULATION} }[$b] -> { LOCATION }[1]; # yes, same method for string access
my $trait = @{ $input ->{POPULATION} }[$b] -> { TRAIT };
$contag_events++;
next;
}
}
}
return $input;
}
#--------------------------------------------------------------------
sub migrate() # rather similar to reproduce
{
my $input = $_[0];
my $migrats = 0; # how many
for(my $b=0;$b<=$popsize-1;$b++)
{
my $rand = rand();
if($rand < $mig_rate) # just migrate every 10th one
{ # choose
$rand = rand(); # random direction
my $xcoord = @{ $input ->{POPULATION} }[$b] ->{LOCATION}[0];
my $ycoord = @{ $input ->{POPULATION} }[$b] ->{LOCATION}[1];
if($rand <= 0.125 && $xcoord >0 && $ycoord >0) # go up and left
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[0]--;
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[1]--;
$pop[$xcoord][$ycoord]--;
$pop[$xcoord-1][$ycoord-1]++;
$migrats++;
}
elsif($rand <= 0.25 && $xcoord >0) # go straight up
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[0]--;
$pop[$xcoord-1][$ycoord]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
elsif($rand <= 0.375 && $xcoord >0 && $ycoord <$cols-1) # go up and right
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[0]--;
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[1]++;
$pop[$xcoord-1][$ycoord+1]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
elsif($rand <= 0.5 && $ycoord <$cols-1) # go straight right
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[1]++;
$pop[$xcoord][$ycoord+1]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
elsif($rand <= 0.625 && $xcoord <$line-1 && $ycoord <$cols-1) # go down and right
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[0]++;
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[1]++;
$pop[$xcoord+1][$ycoord+1]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
elsif($rand <= 0.75 && $xcoord <$line-1) # go straight down
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[0]++;
$pop[$xcoord+1][$ycoord]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
elsif($rand <= 0.825 && $xcoord < $line-1 && $ycoord > 0) # go down and left
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[0]++;
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[1]--;
$pop[$xcoord+1][$ycoord-1]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
elsif($ycoord >0) # go straight left
{
@{ $input ->{POPULATION} }[$b] ->{LOCATION}[1]--;
$pop[$xcoord][$ycoord-1]++;
$pop[$xcoord][$ycoord]--;
$migrats++;
}
}
}
return $input;
}
#--------------------------------------------------------------------------------------------#
sub evenly_initialise() # initialise population
{
my @seq_array; # initialise array for seqs
my $indiv_fitness = 0; # fitness of string, only used when seln. operating
my $av_fitness = 0; # for whole pop, only used when seln. operating
my $st_dev = 0; # ditto
for(my $x=0; $x<=$line-1; $x++)
{
for(my $y=0; $y<=$cols-1; $y++)
{
if($popsize < ($line*$cols)) { die "popsize too small" };
for(my $z=0;$z<=($popsize/($line*$cols))-1;$z++) # run through pop
{
my @one_seq; # re-initialise each string
my $rand = rand(); # decide on horiz trait state
my $horiz_trait;
my $initial_id;
if ($rand < 0.25)
{
$horiz_trait = "A";
}
elsif ($rand < 0.5)
{
$horiz_trait = "B";
}
elsif ($rand < 0.75)
{
$horiz_trait = "C";
}
else
{
$horiz_trait = "O";
}
$rand = rand(); # decide on vertical trait state
if ($rand < 0.25)
{
$initial_id = 1;
}
elsif ($rand < 0.5)
{
$initial_id = 2;
}
elsif ($rand < 0.75)
{
$initial_id = 3;
}
else
{
$initial_id = 4;
}
# RECORD OF SEQUENCE
my $new_seq =
{
LOCATION => [ $x, $y ], # all are initially in square 0,0
SERIAL_NUM => $initial_id, # for lineage mapping
TRAIT => $horiz_trait, # will pass by contagion
};
push @seq_array, $new_seq;
}
}
}
# RECORD OF POP
my $new_pop =
{
POPULATION => [ @seq_array ],
FITNESS => $av_fitness, # if required
STDEV => $st_dev, # if required
};
return $new_pop;
}
#============================================================
sub disp_calc() # calculates dispersals
{
my($gene_state, $trait_state) = @_;
my $total_genes = my $total_traits = my $nbrd_genes = my $nbrd_traits =0;
my @genes = @$gene_state;
my @traits = @$trait_state;
for(my $x=0; $x<=$line-1; $x++)
{
for(my $y=0; $y<=$cols-1; $y++)
{
if($genes[$x][$y] >= 0.5)
{
$total_genes++;
unless($x == 0 || $y == 0 || $x == $line-1 || $y == $cols-1)
{
if($genes[$x+1][$y] >= 0.5
|| $genes[$x-1][$y] >= 0.5
|| $genes[$x][$y+1] >= 0.5
|| $genes[$x][$y-1] >= 0.5
|| $genes[$x+1][$y+1] >= 0.5
||$genes[$x-1][$y-1] >= 0.5
|| $genes[$x+1][$y-1] >= 0.5
|| $genes[$x-1][$y+1] >= 0.5)
{
$nbrd_genes++;
}
}
else
{
if($x == 0 && $y == 0) # top left
{
if($genes[$x+1][$y] >= 0.5
|| $genes[$x][$y+1] >= 0.5
|| $genes[$x+1][$y+1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($x == 0 && $y == $cols-1) # top right
{
if($genes[$x+1][$y] >= 0.5
|| $genes[$x][$y-1] >= 0.5
|| $genes[$x+1][$y-1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($x == $line-1 && $y == 0) # bottom left
{
if($genes[$x-1][$y] >= 0.5
|| $genes[$x][$y+1] >= 0.5
|| $genes[$x-1][$y+1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($x == $line-1 && $y == $cols-1) # bottom right
{
if($genes[$x-1][$y] >= 0.5
|| $genes[$x][$y-1] >= 0.5
|| $genes[$x-1][$y-1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($x == 0) # top row
{
if($genes[$x+1][$y] >= 0.5
|| $genes[$x][$y+1] >= 0.5
|| $genes[$x][$y-1] >= 0.5
|| $genes[$x+1][$y+1] >= 0.5
|| $genes[$x+1][$y-1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($y == 0) # left column
{
if($genes[$x+1][$y] >= 0.5
|| $genes[$x-1][$y] >= 0.5
|| $genes[$x][$y+1] >= 0.5
|| $genes[$x+1][$y+1] >= 0.5
|| $genes[$x-1][$y+1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($x == $line-1) # bottom row
{
if($genes[$x-1][$y] >= 0.5
|| $genes[$x][$y+1] >= 0.5
|| $genes[$x][$y-1] >= 0.5
||$genes[$x-1][$y-1] >= 0.5
|| $genes[$x-1][$y+1] >= 0.5)
{
$nbrd_genes++;
}
}
elsif($y == $cols-1) # right column
{
if($genes[$x+1][$y] >= 0.5
|| $genes[$x-1][$y] >= 0.5
|| $genes[$x][$y-1] >= 0.5
||$genes[$x-1][$y-1] >= 0.5
|| $genes[$x+1][$y-1] >= 0.5)
{
$nbrd_genes++;
}
}
}
}
if($traits[$x][$y] >= 0.5)
{
$total_traits++;
unless($x == 0 || $y == 0 || $x == $line-1 || $y == $cols-1)
{
# print "\n$x, $y not on fringe";
if($traits[$x+1][$y] >= 0.5
|| $traits[$x-1][$y] >= 0.5
|| $traits[$x][$y+1] >= 0.5
|| $traits[$x][$y-1] >= 0.5
|| $traits[$x+1][$y+1] >= 0.5
||$traits[$x-1][$y-1] >= 0.5
|| $traits[$x+1][$y-1] >= 0.5
|| $traits[$x-1][$y+1] >= 0.5)
{
$nbrd_traits++;
}
}
else
{
if($x == 0 && $y == 0) # top left
{
if($traits[$x+1][$y] >= 0.5
|| $traits[$x][$y+1] >= 0.5
|| $traits[$x+1][$y+1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($x == 0 && $y == $cols-1) # top right
{
if($traits[$x+1][$y] >= 0.5
|| $traits[$x][$y-1] >= 0.5
|| $traits[$x+1][$y-1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($x == $line-1 && $y == 0) # bottom left
{
if($traits[$x-1][$y] >= 0.5
|| $traits[$x][$y+1] >= 0.5
|| $traits[$x-1][$y+1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($x == $line-1 && $y == $cols-1) # bottom right
{
if($traits[$x-1][$y] >= 0.5
|| $traits[$x][$y-1] >= 0.5
|| $traits[$x-1][$y-1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($x == 0) # top row
{
if($traits[$x+1][$y] >= 0.5
|| $traits[$x][$y+1] >= 0.5
|| $traits[$x][$y-1] >= 0.5
|| $traits[$x+1][$y+1] >= 0.5
|| $traits[$x+1][$y-1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($y == 0) # left column
{
if($traits[$x+1][$y] >= 0.5
|| $traits[$x-1][$y] >= 0.5
|| $traits[$x][$y+1] >= 0.5
|| $traits[$x+1][$y+1] >= 0.5
|| $traits[$x-1][$y+1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($x == $line-1) # bottom row
{
if($traits[$x-1][$y] >= 0.5
|| $traits[$x][$y+1] >= 0.5
|| $traits[$x][$y-1] >= 0.5
||$traits[$x-1][$y-1] >= 0.5
|| $traits[$x-1][$y+1] >= 0.5)
{
$nbrd_traits++;
}
}
elsif($y == $cols-1) # right column
{
if($traits[$x+1][$y] >= 0.5
|| $traits[$x-1][$y] >= 0.5
|| $traits[$x][$y-1] >= 0.5
||$traits[$x-1][$y-1] >= 0.5
|| $traits[$x+1][$y-1] >= 0.5)
{
$nbrd_traits++;
}
}
}
}
}
}
my $g_dispersal;
my $t_dispersal;
return ($total_genes-$nbrd_genes, $total_traits-$nbrd_traits);
}
© Copyright Journal of Artificial Societies and Social Simulation, [2002]