program Culture;
{Robert Axelrod, School of Public Policy, U Michigan, Ann Arbor, MI 48109}
{ This version is based on version 2.6 Axelrod's Org program for his Cultural Model.}
const
Version = 2.6; {Program Version Number}
test = true; {if test=true, no run_number is read or written.}
{Run # is set to 0 for test=true runs}
old_random_seed = 0; {0 means new seed, else enter an old seed to reuse}
cyclemax_float = 10; {# cycles in each reporting period. Cylce is one active actor, ie event. eg 1000}
periodmax = 10; {number of periods of cyclemax each.Max 30,000 }
mapmax = 1; {number of maps}
popmax = 1; {number of replications of each map, each with cyclemax*periodmax active actors}
mutuation_rate_perpop_percycle = 0.0; {mutation rate, per pop, per cycle for drift, eg 0.0}
broadcast_rate_per_10000_cycles = 0; {cylces in 10,000 that are interactions with broadcaster, not neighbors}
{ e.g. 200 means 2% of cycles are from broadcast, 98% from a neighbor}
flat_map = true; {true is old method of no wrap around. false is torus}
{production parameters}
N_Xmax = 1; {Map sizes used. (max 5. See Initialize_Run, and def production_array)}
Xmax_tops = 50; { largest allowable value of Xmax1, Xmax2, etc.}
Ymax_tops = Xmax_tops; {assume square map}
Xmax1 = 10; { first map size East to West: Square map assumed. Max is Xmax_tops}
Xmax2 = 5; { second map size, etc.}
Xmax3 = 7;
Xmax4 = 12;
Xmax5 = 17;
N_interaction_types = 1; {neighborhood types used}
Interaction_types_tops = 16;{ largerst allowable is 16}
Interaction_Types1 = 4; {first used. : 4=NEWS, 8=3x3, 12=diamond, 16=diamond+NEWS repeated}
Interaction_Types2 = 8;
Interaction_Types3 = 12;
Interaction_Types4 = 0;
Interaction_Types5 = 0;
N_bitmax = 1; {Chromosome sizes used (max 5)}
bitmax_tops = 15; { largest allowable bitmax1, bitmax2, etc}
bitmax1 = 5; { Size of chromosome in first set of pops, ie bits in indiv's culture. Max bitmax_tops}
bitmax2 = 10; { second string length, etc.}
bitmax3 = 15;
bitmax4 = 0;
bitmax5 = 0;
N_allelemax = 1; {Alleles in indiv's culture (number of different values to be used) (max 5)}
allelemax1 = 15; {Alleles in each gene in first set of pops}
allelemax2 = 10; {Alleles in each gene in first set of pops}
allelemax3 = 15;
allelemax4 = 0;
allelemax5 = 0;
{reporting options: }
TextGridReport = false; {Report period grid data to text window; slows things down}
Graphics = false; {Show graphics in drawing window if true}
Graph_distance = false; {Plot cul. distance in drawing window if true (and graphics=true)}
WritePeriodicStats = true; {Writes period statistics to period_output if true}
WriteRegionStats = false; {Writes data on regions each period, huge unless few periods}
WriteRegionStatsAtEnd = true; {Writes data on regions at end of a case to Region_Output if neighbor_clear}
stack_map = true; {if false, can see multiple maps in drawing window}
report_bit_origin = false; {write origin of first bit as Xmax*(y-1)+x-1}
{graphics view: }
w = 10; {Graphics constant: width of a cell, eg 20}
h_offset = 20; {Graphics constant: where grid starts in drawing window, horizontal offset, eg 20}
v_offset = 20; {Graphics constant: where grid starts, eg 20}
pen_size = 4; {Graphics constant: height and width of pen, eg 6}
{compiler constants}
Xmax_topsPlus1 = Xmax_tops + 1; {highest values allowed of Xmax}
Ymax_topsPlus1 = Ymax_tops + 1;
XmaxYmax_tops = Xmax_tops * Ymax_tops;
Run_Number_File_Name = 'Run Number File';
type
x_type = 1..xmax_tops; {range of indiv's x coord}
y_type = 1..ymax_tops; {range of indiv's y coord}
bit_type = 1..bitmax_tops; {bit on culture}
culture_array = array[0..Xmax_topsPlus1, 0..Ymax_topsPlus1, 1..bitmax_tops] of integer; {for x, y,bit: max sizes allowed}
culture_array_pointer = ^culture_array;
Culture_Array_Handle = ^culture_array_pointer;
done_array = array[0..Xmax_topsPlus1, 0..Ymax_topsPlus1] of boolean;{for analyzing regions. max actual geo }
done_array_pointer = ^done_array;
done_array_handle = ^done_array_pointer;
direction_type = 1..16; {neighbors: see pop initialization. Starts with north, e, w, s order}
direction_array = array[1..16] of integer; {for x and y moves}
distance_array = array[1..Xmax_tops, 1..Ymax_tops] of integer; {for distances right and down}
distance_array_pointer = ^distance_array;
distance_array_handle = ^distance_array_pointer;
production_array = array[1..5] of integer; {for storing production values}
location_array = array[1..XmaxYmax_tops] of integer; {for x and y loc of first member of a region}
location_array_pointer = ^location_array;
location_array_handle = ^location_array_pointer;
var
initial_datetime, end_datetime: datetimerec; {date, etc.}
initial_hour, end_hour: longint;
start_time, end_time, duration: longint; {for calc of run's duration}
random_seed: integer;
run_number: integer;
datafile: text; {input data}
final_output: text;
region_output: text;
period_output: text;
pop: integer; {current population number, 1...popmax}
period: longint; {current period number, 1..periodmax}
ix, iy: integer; {x,y coordinate of current active individual}
jx: integer; {x coord of a neighbor, 0..xmax+1}
jy: integer; {y coord of a neighbor, 0..ymax+1}
culture: Culture_Array_Handle; {culture of ix,iy,bit - an integer variable}
{-1 for beyond bord, or 1... allelemax}
{x=0 is beyond internal left border, x=xmax+1 is beyond right border, etc.}
{This allows indivs on internal border to get automatic non-matches with illegal neighbors}
xmove, ymove: direction_array; {jointly defines North, East, West, South direction}
bit: integer; {the current bit of the culture}
neighbor: direction_type; {neighbor, 1 to 16}
direction: direction_type; {16 neighbors; starts with north, s, e, w order, then NE etc.}
NoMatchPeriodCount: integer; {count of no matches in this period}
NoMatchPopCount: integer; {count of no matches in this pop}
match: boolean; {true if active matches selected neighbor on selected bit}
bit_distance_sum: longint; {sum (in current cycle) of neighbors' bit distances over pop in a pop, }
n_zero_distance: integer; {num (in current cycle) of boundaries of zero bit distance}
talk_count: integer; {counts number of interactions tried in this cycle}
r_dis: distance_array_handle; {bit distance of an actor to its right neighbor}
d_dis: distance_array_handle; {bit distance of an actor to its down neighbor}
iq: integer; {index of ready queue in analyze regions}
done: done_array_handle; {true if the actor has been processed by regional analysis}
regions: integer; {count of number of regions}
region_size: integer; {number actors in current region}
x_region, y_region: location_array_handle; {x and y loc of first member of a region }
Xmax: 1..xmax_tops; {size of map east to west}
Ymax: 1..ymax_tops; {size of map north to south}
interaction_types: 1..Interaction_types_tops; {neighborhood type: 4=NEWS, 8=3x3, 12=diamond, 16=diamond+NEWS repeated}
bitmax: 1..bitmax_tops; {number of genes in culture string}
allelemax: integer; {number of alleles in each gene}
Xmax_value: production_array; {Xmax's first, second, third, etc values}
Interaction_Types_value: production_array; {interaction_type first, second, third etc values}
bitmax_value: production_array; {bitmax's first, second, third, etc values}
allelemax_value: production_array; {allelemax's first, second, third, etc values}
XMaxPlus1: integer; {current value of xmax+1}
YMaxPlus1: integer; {current value of ymax+1}
Xmax_index: x_type; {in main program's loop}
bitmax_index: bit_type;
allelemax_index: integer;
XmaxYmax: integer; {current value}
interaction_types_index: integer;
case_number: integer; {sequential count of pops}
broadcast_allowed: boolean; {true if broadcast_rate_per_10000_cycles > 0}
broadcast_cutoff_number: integer; {to get correct prob that an interaction is from a broadcast}
broadcast_count: longint; {number of braodcasts}
Last_X_Boundary: x_type; {for distances. =xmax-1 if flat_world=true, else =xmax}
Last_Y_Boundary: x_type; {for distances. =ymax-1 if flat_world=true, else =ymax}
map: integer; {2=second map used in this set, etc., goes 1 to mapmax}
old_seed: integer; {for repitition of maps}
map_seed: array[1..mapmax] of integer; {seeds to start maps}
interaction_seed: array[1..popmax] of integer; {seeds to start interactions}
map_x_offset: integer; {to displace maps horiztonally, one row for each initial map}
map_y_offset: integer; {to displace maps vertically, one col for each pop type}
origin_first_bit: distance_array_handle; {original site of first bit}
Changes_This_Period: longint; {count of actual number of cultural changes this period}
Zones: integer; {count of number of zones, ie domains separated by maximal cul boundary}
cyclemax: longint; {number of cycles in a period, from cyclemax_float input constant}
{wp: WindowPtr; }
{pointer to my window for drawing, from Org25, discared in Org26}
{ --------------------------------------------------------------- }
function random_one_to_n (n: longint): longint;
{proc returns a random number between 1 and n; modified from Bennett's random_range}
var
ub, lb: integer;
r: integer;
begin {random gives # betw -32768 and 32767}
ub := 32767 - (32767 mod n);
lb := -32768 - (-32768 mod n); {truncate distrib on 2 ends so that later mod is OK}
repeat
r := random;
until (r <= ub) and (r >= lb); {make sure random genrated is in truncated (even) distrib}
random_one_to_n := abs(r mod n) + 1;
end; {random function}
{ --------------------------------------------------------------- }
function Random0or1: integer;
{ returns 0 or 1 with equal probability}
begin
if random < 0 then {random gives # betw -32768 and 32767}
Random0or1 := 0
else
Random0or1 := 1;
end; {Random0or1 fcn}
{ --------------------------------------------------------------- }
function ipoisn (mean: real): integer;
{returns a random integer from a Poisson distribution.}
{The mean of the distribution is passed into the function. }
{ Adapted from S. Forrest, GA Program v3., 6/6/83, p8}
{ Lamda not inputted since always calculated in this fcn}
var
k: integer; {k term in Poisson disribution}
cum, {current value of probability density fcn}
term, {current value of prob distibu? fcn}
lambda_term, {cacl to be exp(-mean)}
tempnum: real; {random number beteween 0 and 1 from uniform distibution}
begin
if mean < 0 then
write(' Error from ipoisn fcn, mean must be non-negative')
else
begin
lambda_term := exp(-1 * mean);
tempnum := (random + 32768) / 65635; {random gives # betw -32768 and 32767}
{tempnum gives uniform between 0 and 1 }
{Here, the cum density fcn for successive values of k is computed. When the cdf}
{becomes greater than the random variable, tempnum, we quit. Since k must be}
{non-negative and since prob(k=0) = lambda_term, k will start out at 0 with}
{the other values initialized acordingly.}
k := 0;
term := lambda_term;
cum := lambda_term;
while (cum <= tempnum) and (term > 0.0000005) do {cutoff for high end of distribution}
begin
k := succ(k);
term := (term * mean) / k;
cum := cum + term;
end;{cylce of cumulation}
ipoisn := k; {return k as the random number}
end; {else}
{writeln('ipoisn test: mean, output k ', mean : 8 : 3, k : 3);}
end; {ipoisn fcn}
{ --------------------------------------------------------------- }
procedure create_graphics_window;
const
size = 10; {font}
var
style: integer;
wrect: rect;
begin
wrect.left := 0;
wrect.top := 40;
wrect.right := 500;
if h_offset + xmax_tops * w > 450 then
wrect.right := h_offset + xmax_tops * w + 50;
if Graph_distance then
wrect.bottom := v_offset + (ymax_tops * w + 10) * popmax + 150 + periodmax
else
wrect.bottom := v_offset + (ymax_tops * w + 10) * popmax + 100;
setdrawingrect(wrect); {discarded in Org25, restored in Org26}
showdrawing;
{following from Org25, dicarded inOrg26}
{wp := NewWindow(nil, wrect, 'Map Window', true, 0, WindowPtr(-1), true, 0);}
{setport(wp);}
TextSize(size);
getfnum('Courier', style);
textfont(style);
end;
{ --------------------------------------------------------------- }
procedure initialize_graphics_window;
begin
moveto(h_offset, 10);
if test then
drawstring(concat(' Ax Org Program, version ', stringof(version : 3 : 2)))
else
drawstring(concat(' Ax Org Program, version ', stringof(version : 3 : 2), ', Run', stringof(run_number : 4), '.'));
moveto(h_offset, 20);
drawstring(concat(' This run begun on ', stringof(initial_datetime.month : 2), '/', stringof(initial_datetime.day : 2), '/', stringof(initial_datetime.year : 4), ' at ', stringof(initial_datetime.hour : 2), ':', stringof(initial_datetime.minute : 2), '.'));
end;
{ --------------------------------------------------------------- }
procedure frame_rectangle; {for outline of map}
begin
pensize(1, 1);
penNormal;
framerect(v_offset + 29 + pen_size + map_y_offset, h_offset - 1 + map_x_offset, ymax * w + v_offset + 31 + map_y_offset, xmax * w - pen_size + h_offset + 1 + map_x_offset);
pensize(pen_size, pen_size);
end;
{ --------------------------------------------------------------- }
procedure write_run_info;
begin
rewrite(final_output, 'Final_Output'); {open output of final file for writing to Excel}
rewrite(region_output, 'Region_Output'); {open output of region data for writing to Excel}
if WritePeriodicStats then
rewrite(period_output, 'Period_Output '); {open output of period data}
gettime(initial_datetime);
initial_hour := initial_datetime.hour; {to force long int}
start_time := 60 * 60 * initial_hour + 60 * initial_datetime.minute + initial_datetime.second;
Writeln('Output of Axelrod''s Org Program, Version ', Version : 5 : 2);
Write(' This run begun on ', initial_datetime.month : 2, '/', initial_datetime.day : 2);
Writeln('/', initial_datetime.year : 4, ' at ', initial_datetime.hour : 2, ':', initial_datetime.minute : 2, '.');
Writeln(cyclemax : 5, ' cycles in each reporting period');
Writeln(periodmax : 5, ' maximum number of periods in each pop');
Writeln(mapmax : 5, ' number of maps');
Writeln(popmax : 5, ' number of pops for each map');
Writeln(broadcast_rate_per_10000_cycles : 5, ' broadcast rate per 10,000 cycles');
Writeln(flat_map : 5, ' flat map. True means no wrap around, false means torus map');
{ writeln(Interaction_Types : 5, ' interaction types . 4 means 4 neighbors , etc . See documentation . ');}
Writeln(mutuation_rate_perpop_percycle : 7 : 3, ' mutation rate, per pop, per cycle');
Writeln(final_output, 'Output of Axelrod''s Org Program, Version ', Version : 5 : 2);
Write(final_output, ' This run begun on ', initial_datetime.month : 2, '/', initial_datetime.day : 2);
Writeln(final_output, '/', initial_datetime.year : 4, ' at ', initial_datetime.hour : 2, ':', initial_datetime.minute : 2, '.');
Writeln(final_output, cyclemax : 5, ' cycles in each reporting period');
Writeln(final_output, periodmax : 5, ' maximum number of periods in each pop');
Writeln(final_output, mapmax : 5, ' number of maps');
Writeln(final_output, popmax : 5, ' number of pops for each map');
Writeln(final_output, broadcast_rate_per_10000_cycles : 5, ' broadcast rate per 10,000 cycles');
Writeln(final_output, flat_map : 5, ' flat map. True means no wrap around, false means torus map');
{ writeln(Interaction_Types : 5, ' interaction types. 4 means 4 neighbors, etc. See documentation.');}
Writeln(final_output, mutuation_rate_perpop_percycle : 7 : 3, ' mutation rate, per pop, per cycle');
Writeln(region_output, 'Region Output of Axelrod''s Org Program, Version ', Version : 5 : 2);
Write(region_output, ' This run begun on ', initial_datetime.month : 2, '/', initial_datetime.day : 2);
Writeln(region_output, '/', initial_datetime.year : 4, ' at ', initial_datetime.hour : 2, ':', initial_datetime.minute : 2, '.');
Writeln(period_output, 'Period Output of Axelrod''s Org Program, Version ', Version : 5 : 2);
Write(period_output, ' This run begun on ', initial_datetime.month : 2, '/', initial_datetime.day : 2);
Writeln(period_output, '/', initial_datetime.year : 4, ' at ', initial_datetime.hour : 2, ':', initial_datetime.minute : 2, '.');
if old_random_seed = 0 then
begin {generate new seed}
random_seed := initial_datetime.hour + initial_datetime.minute + initial_datetime.second + (initial_datetime.second * 300) + (initial_datetime.minute * initial_datetime.hour) + (initial_datetime.minute * initial_datetime.second);
randseed := random_seed;
Writeln(' New random seed ', randseed : 6, '.');
Writeln(final_output, ' Final Output.', ' New random seed', randseed : 8, '.');
Writeln(region_output, ' New random seed', randseed : 8, '.');
if WritePeriodicStats then
Writeln(period_output, ' New random seed', randseed : 8, '.');
end
else
begin {use old seed, which was inputed as constant}
randseed := old_random_seed;
Writeln(' Old random seed ', randseed : 6, '.');
Writeln(final_output, ' Old random seed', randseed : 8, '.');
Writeln(region_output, ' Old random seed', randseed : 8, '.');
if WritePeriodicStats then
Writeln(period_output, ' Old random seed', randseed : 8, '.');
end;
if test then {run number}
run_number := 0
else
begin {not a test; read run #}
reset(datafile, Run_Number_File_Name);
readln(datafile, run_number);
run_number := run_number + 1;
close(datafile);
rewrite(datafile, Run_Number_File_Name);
writeln(datafile, run_number);
writeln('Run number', run_number : 4, '.');
close(datafile);
Writeln(final_output, 'Run number', run_number : 4, '.');
Writeln(region_output, 'Run number', run_number : 4, '.');
if WritePeriodicStats then
Writeln(period_output, 'Run number', run_number : 4, '.');
end; {if for run number }
end;
{ --------------------------------------------------------------- }
procedure set_up_production_arrays;
var
itemp: integer;
begin
Xmax_value[1] := Xmax1; {assume square map, so Ymax set in main program}
Xmax_value[2] := Xmax2;
Xmax_value[3] := Xmax3;
Xmax_value[4] := Xmax4;
Xmax_value[5] := Xmax5;
for itemp := 1 to N_xmax do
begin {check enough space reserved for the map}
if Xmax_tops < Xmax_value[itemp] then
begin
writeln('**** Halted because Xmax_tops not set big enough ****');
writeln(final_output, ' * * * * Halted because Xmax_tops not set big enough * * * * ');
moveto(h_offset, 30);
drawstring(' * * * * Halted because Xmax_tops not set big enough * * * * ');
halt;
end;
end;
Interaction_Types_value[1] := Interaction_types1;
Interaction_Types_value[2] := Interaction_types2;
Interaction_Types_value[3] := Interaction_types3;
Interaction_Types_value[4] := Interaction_types4;
Interaction_Types_value[5] := Interaction_types5;
bitmax_value[1] := bitmax1;
bitmax_value[2] := bitmax2;
bitmax_value[3] := bitmax3;
bitmax_value[4] := bitmax4;
bitmax_value[5] := bitmax5;
allelemax_value[1] := allelemax1;
allelemax_value[2] := allelemax2;
allelemax_value[3] := allelemax3;
allelemax_value[4] := allelemax4;
allelemax_value[5] := allelemax5;
end;
{ --------------------------------------------------------------- }
procedure set_up_move_definitions;
begin
Xmove[1] := 0; {define North, in seeking neighbors}
Ymove[1] := -1;
Xmove[2] := 1; {define East}
Ymove[2] := 0;
Xmove[3] := -1; {define West}
Ymove[3] := 0;
Xmove[4] := 0; {define South}
Ymove[4] := 1;
Xmove[5] := 1; {define NE}
Ymove[5] := -1;
Xmove[6] := 1; {define SE}
Ymove[6] := 1;
Xmove[7] := -1; {define SW}
Ymove[7] := 1;
Xmove[8] := -1; {define NW}
Ymove[8] := -1;
Xmove[1 + 8] := 0; {define 2 North}
Ymove[1 + 8] := -2;
Xmove[2 + 8] := 2; {define 2 East}
Ymove[2 + 8] := 0;
Xmove[3 + 8] := -2; {define 2 West}
Ymove[3 + 8] := 0;
Xmove[4 + 8] := 0; {define 2 South}
Ymove[4 + 8] := 2;
Xmove[1 + 12] := 0; {repeat for directions 13-16: repeat North}
Ymove[1 + 12] := -1;
Xmove[2 + 12] := 1; {repeat East}
Ymove[2 + 12] := 0;
Xmove[3 + 12] := -1; {repeat West}
Ymove[3 + 12] := 0;
Xmove[4 + 12] := 0; {repeat South}
Ymove[4 + 12] := 1;
end;
{ --------------------------------------------------------------- }
procedure Initialize_run;
var
temp: integer; {for tab index}
begin
cyclemax := round(cyclemax_float);
culture := Culture_Array_Handle(NewHandle(sizeof(Culture_Array)));
done := done_Array_Handle(NewHandle(sizeof(done_Array)));
r_dis := distance_Array_Handle(NewHandle(sizeof(distance_Array)));
d_dis := distance_Array_Handle(NewHandle(sizeof(distance_Array)));
X_region := location_Array_Handle(NewHandle(sizeof(location_Array)));
Y_region := location_Array_Handle(NewHandle(sizeof(location_Array)));
origin_first_bit := distance_Array_Handle(NewHandle(sizeof(distance_Array)));
if broadcast_rate_per_10000_cycles > 0 then
begin
broadcast_count := 0;
broadcast_allowed := true; {next line: to compare to -32768 to 32767}
broadcast_cutoff_number := round((broadcast_rate_per_10000_cycles / 10000.0) * 65536.0 - 32768.0);
{writeln('broadcast_cutoff_number', broadcast_cutoff_number : 6); }
end;
write_run_info; {write run parameters to files}
write('Case ');
if mapmax <> 1 then {repeating the same map}
Write('Map ');
writeln('Pop Xmax Neigh String Alleles Clear Period Xmax Regions tot dis # zero dis');
write(final_output, 'Case ');
if mapmax <> 1 then
Write(final_output, 'Map ');
writeln(final_output, 'Pop Xmax Neigh String Alleles Clear Period Xmax Regions tot dis # zero dis');
if WriteRegionStatsAtEnd = true then
writeln(region_output, 'Case Map Pop Region X Y Size Cul....then Distance to prev region');
if graphics then
begin
create_graphics_window;
initialize_graphics_window
end;
if WritePeriodicStats then
begin
writeln(period_output, 'case map pop period cycle dist #0 dist regions changes zones');
end;
set_up_production_arrays;
if Interaction_Types_Tops > 16 then
begin
writeln('*** Error in setting contant: Interaction_Types_Tops > 16. Program halted. ***');
halt;
end;
set_up_move_definitions;
end; {initialize run procedure}
{ --------------------------------------------------------------- }
procedure write_graphics_report;
begin
frame_rectangle; {to restore outer edge}
forecolor(whitecolor); {update text at bottom}
paintrect(10 + v_offset, h_offset, 30 + v_offset, 180 + h_offset);
forecolor(blackcolor);
moveto(h_offset, 18 + v_offset);
drawstring(concat('Period =', stringof(period : 5)));
if cyclemax <> 1000 then
drawstring(concat(' (Each period is', stringof(cyclemax : 8), ' cycles.)'));
if Graph_distance then
begin
moveto(h_offset, w * ymax_tops + 30 + v_offset);
drawstring(concat('Sum of distances =', stringof(bit_distance_sum : 4)));
moveto(h_offset, w * ymax_tops + 50 + v_offset + period);
penNormal;
line(bit_distance_sum, 0); {draw horizontal line}
pensize(pen_size, pen_size); {restore size of pen}
end; {graph cultural distance}
end;
{ --------------------------------------------------------------- }
procedure Initialize_pop;
var
x, y: integer; {local variable for coord of an individual}
bit: integer; {local variable for number of gene on the cultrual chormosome}
begin
if TextGridReport then
writeln('Pop ', pop : 3, ':');
NoMatchPopCount := 0; {count of times no match is found}
for x := 1 to xmax do
begin
for y := 1 to ymax do
begin
for bit := 1 to bitmax do
begin {initial culture for internal places}
culture^^[x, y, bit] := Random_one_to_n(allelemax) - 1; {0..allelemax-1, sarting org13}
end;
origin_first_bit^^[x, y] := Xmax * (y - 1) + x - 1; {record origin of the first bit, 00 to 99 for 10x10 map}
end;
end; {x loop}
if stack_map then
begin
map_x_offset := 0; {draw maps on top of each other}
map_y_offset := 0
end
else
begin {write maps in rows for same initial conditions}
map_x_offset := (w * xmax + 10) * (pop - 1);
map_y_offset := (w * ymax + 10) * (map - 1);
end;
end; {initialize_pop procedure}
{ --------------------------------------------------------------- }
procedure Initialize_period;
begin
Changes_This_Period := 0; {count of number of actual cultural changes in this period}
end; {initialize_period procedure}
{ --------------------------------------------------------------- }
procedure report_change_test1; {from interact procedure. reports first 5 bits.}
begin
writeln('report_converge_test1.match=', match : 2, ' ix=', ix : 3, '. iy=', iy : 3, ' ini cul = ', culture^^[ix, iy, 1] : 1, culture^^[ix, iy, 2] : 1, culture^^[ix, iy, 3] : 1, culture^^[ix, iy, 4] : 1, culture^^[ix, iy, 5] : 1);
{writeln('report_converge_test1. jx=', jx : 3, '. jy=', jy : 3, ' init cul = ', culture^^[jx, jy, 1] : 1, culture^^[jx, jy, 2] : 1, culture^^[jx, jy, 3] : 1, culture^^[jx, jy, 4] : 1, culture^^[jx, jy, 5] : 1);}
end;
{ --------------------------------------------------------------- }
procedure report_change_test2; {from interact procedure}
begin
writeln('report_change_test2. ix=', ix : 3, '. iy=', iy : 3, 'js new bit=', bit : 2, '. new cul = ', culture^^[ix, iy, 1] : 1, culture^^[ix, iy, 2] : 1, culture^^[ix, iy, 3] : 1, culture^^[ix, iy, 4] : 1, culture^^[ix, iy, 5] : 1);
{ writeln('report_change_test2. jx=', jx : 3, '. jy=', jy : 3, 'js new bit=', bit : 2, '. new cul = ', culture^^[jx, jy, 1] : 1, culture^^[jx, jy, 2] : 1, culture^^[jx, jy, 3] : 1, culture^^[jx, jy, 4] : 1, culture^^[jx, jy, 5] : 1);}
end;
{ --------------------------------------------------------------- }
procedure change_if_possible_toward_neigh; {i converges one bit to j if possible}
var
try_another_bit: boolean; {control on whether still searching for another bit that differs}
bit_count: integer; {count so as to give up when all bit tried, ie x=y}
bit_try: integer; {bit being tried looking for dissimilarity}
begin
{ report_change_test1;}
bit_try := random_one_to_n(bitmax); {initial bit location tried}
try_another_bit := true;
bit_count := 1;
repeat
if culture^^[ix, iy, bit_try] <> culture^^[jx, jy, bit_try] then
begin
culture^^[ix, iy, bit_try] := culture^^[jx, jy, bit_try]; {i converges since unequal on current bit}
try_another_bit := false; { done with search }
Changes_This_Period := Changes_This_Period + 1; {to count actual number of changes, cf interact from broadcast}
if bit_try = 1 then {update where the first bit came from}
origin_first_bit^^[ix, iy] := origin_first_bit^^[jx, jy];
end
else
begin
bit_count := bit_count + 1;
bit_try := (bit_try + 1) mod bitmax + 1;
if bit_count > bitmax then {give up because just tried all bits}
try_another_bit := false;
end;
until try_another_bit = false;
{ report_change_test2;}
end;{interact procedure}
{------------------------------------------------------------}
procedure interact_from_neighbor;
begin
repeat {get an internal direction}
direction := random_one_to_n(Interaction_Types); {select location for interaction}
if flat_map then
begin
jx := ix + xmove[direction]; {cacl coords of selected neighbor, j}
jy := iy + ymove[direction]
end
else {torus map}
begin {note: have to add Xmax or Ymax to be sure to get positive number}
jx := (ix + xmove[direction] - 1 + Xmax) mod xmax + 1; {cacl coords of selected neighbor, j}
jy := (iy + ymove[direction] - 1 + Ymax) mod ymax + 1;
end;
until ((flat_map = false) or ((jx > 0) and (jx < XMaxPlus1) and (jy > 0) and (jy < YMaxPlus1)));
{other is in legal area}
match := (culture^^[ix, iy, bit] = culture^^[jx, jy, bit]); {match is true or false}
{writeln('test3. match=', match : 3, ' iy=', iy : 4, '. y direction = ', direction : 4, ' . ymove = ', ymove[direction] : 4);}
if match then {if the chosen bit of the selected neighbor matches then interact }
change_if_possible_toward_neigh {ie do only one interaction, as in Org 06}
end;
{ --------------------------------------------------------------- }
procedure change_if_possible_toward_broadcast; {i converges one bit to j if possible}
var
try_another_bit: boolean; {control on whether still searching for another bit that differs}
bit_count: integer; {count so as to give up when all bit tried, ie x=y}
bit_try: integer; {bit being tried looking for dissimilarity}
begin
{report_change_test1;}
bit_try := random_one_to_n(bitmax); {initial bit location tried}
try_another_bit := true;
bit_count := 1;
repeat
if culture^^[ix, iy, bit_try] <> 0 then
begin
culture^^[ix, iy, bit_try] := 0; {i converges since unequal on current bit}
Changes_This_Period := Changes_This_Period + 1; {to count actual number of changes, cf change to neighbor}
try_another_bit := false; { done with search }
end
else
begin
bit_count := bit_count + 1;
bit_try := (bit_try + 1) mod bitmax + 1;
if bit_count > bitmax then {give up because just tried all bits}
try_another_bit := false;
end;
until try_another_bit = false;
{report_change_test2;}
end;{interact procedure}
{------------------------------------------------------------}
procedure interact_from_broadcast;
begin
{broadcast_count := broadcast_count + 1; }
match := culture^^[ix, iy, bit] = 0; {match is true or false}
{writeln('test3. match=', match : 3, ' iy=', iy : 4, '. y direction = ', direction : 4, ' . ymove = ', ymove[direction] : 4);}
if match then {if the chosen bit of the selected neighbor matches then interact }
change_if_possible_toward_broadcast {ie do only one interaction, as in Org 06}
end;
{ --------------------------------------------------------------- }
procedure mutate; {mutate for drift}
var
num_mutations: integer; {number of mutations to perform this cycle}
mut: integer; {index of mutations}
old_allele, new_allele: integer;
begin
num_mutations := ipoisn(mutuation_rate_perpop_percycle); {calc number of mutations}
{writeln(num_mutations : 3, ' n muts'); }
{test}
for mut := 1 to num_mutations do {each mutation}
begin
ix := random_one_to_n(xmax);
iy := random_one_to_n(ymax);
bit := random_one_to_n(bitmax);
Old_allele := culture^^[ix, iy, bit]; {mutate to a DIFFERENT value 1...allelemax}
New_allele := random_one_to_n(allelemax - 1) - 1; {allele-1 corrected ver 1.5; should have been ver 1.3}
if new_allele = old_allele then
new_allele := allelemax; { if initially no change, use this for maxallele}
culture^^[ix, iy, bit] := new_allele;
end;{mutations loop}
end;{mutate procedure}
{ --------------------------------------------------------------- }
procedure report_output_period_test; {from output_period procedure}
var {reports first five bits of culture}
xtemp, ytemp: integer;
begin
for ytemp := 1 to ymax do
begin
write(' y=', ytemp : 3, '. Cul(1..5)=');
for xtemp := 1 to xmax do
begin
write(culture^^[xtemp, ytemp, 1] : 3, culture^^[xtemp, ytemp, 2] : 1, culture^^[xtemp, ytemp, 3] : 1, culture^^[xtemp, ytemp, 4] : 1, culture^^[xtemp, ytemp, 5] : 1);
end;{ixtemp}
writeln;
end; {iytemp}
end;
{ --------------------------------------------------------------- }
procedure set_pen_pat (bit_distance: integer);
begin
if bit_distance = 0 then
penpat(white)
else if bit_distance = 1 then
penpat(ltgray)
else if bit_distance = 2 then
penpat(gray)
else if bit_distance = 3 then
penpat(dkgray)
else
penpat(black);
end;
{ --------------------------------------------------------------- }
procedure draw_right_boundary (xtemp: integer; ytemp: integer);
var
bit_distance: integer;
begin
bit_distance := r_dis^^[xtemp, ytemp];
set_pen_pat(bit_distance);
move(0, pen_size);
line(0, w - 2 * pen_size);
move(w, -w + pen_size);
end;
{ --------------------------------------------------------------- }
procedure draw_bottom_boundary (xtemp: integer; ytemp: integer);
var
bit_distance: integer;
begin
bit_distance := d_dis^^[xtemp, ytemp];
set_pen_pat(bit_distance);
line(w - 2 * pen_size, 0);
if xtemp <= Last_X_Boundary then {draw intersection}
begin {find max distance around interection}
if (r_dis^^[xtemp, ytemp] > bit_distance) then
bit_distance := r_dis^^[xtemp, ytemp];
if (d_dis^^[xtemp + 1, ytemp] > bit_distance) then
bit_distance := d_dis^^[xtemp + 1, ytemp];
if (r_dis^^[xtemp, ytemp + 1] > bit_distance) then
bit_distance := r_dis^^[xtemp, ytemp + 1];
set_pen_pat(bit_distance);
move(pen_size, 0);
line(0, 0); {to draw intersection}
move(pen_size, 0);
end;
end;
{ --------------------------------------------------------------- }
procedure draw_distance_graph;
var
xtemp, ytemp: integer; {geo loc}
begin
moveto(h_offset, 10); {10 is space needed for header on drawing window}
for ytemp := 1 to ymax do
begin
moveto(w + h_offset - pen_size + map_x_offset, w * (ytemp - 1) + v_offset + map_y_offset + 30);
for xtemp := 1 to Last_X_Boundary do
begin
draw_right_boundary(xtemp, ytemp);
end;{xtemp}
moveto(h_offset + map_x_offset, (w) * ytemp + v_offset + map_y_offset + 30);
if (ytemp <= Last_Y_Boundary) then
begin
for xtemp := 1 to xmax do {write row for distances down, only if not last row}
draw_bottom_boundary(xtemp, ytemp); {including intersection}
end; {if not last row}
end; {ytemp}
write_graphics_report;
end;
{ --------------------------------------------------------------- }
procedure CalcDistance; {cultural distances with output}
type
row_array = array[1..xmax_tops] of integer; {for output of distances}
var
xtemp, ytemp: integer; {geo loc}
itemp: integer; {bit position}
X_distance: row_array; {cultural distance between (x,y) to (x+1,y), ie to right }
Y_distance: row_array; {cultural distance between (x,y) to (x,y+1), ie down }
begin
bit_distance_sum := 0; {initialize sum of bit distances to neighbors (over current cycle)}
n_zero_distance := 0; {initialize num of boundaries with zero bit distance}
for ytemp := 1 to ymax do
begin
for xtemp := 1 to xmax do
begin
X_distance[xtemp] := 0; {calcs will be for one row at a time}
Y_distance[xtemp] := 0;
for itemp := 1 to bitmax do
begin {increment distance on x axis, then y axis}
if flat_map then {flat map}
begin
if xtemp <= Last_X_Boundary then {to avoid going off right side, unless torus}
if culture^^[xtemp, ytemp, itemp] <> culture^^[xtemp + 1, ytemp, itemp] then
X_distance[xtemp] := X_distance[xtemp] + 1;
if (ytemp <= Last_Y_Boundary) then {only do down calcs if not last row}
begin
if (culture^^[xtemp, ytemp, itemp] <> culture^^[xtemp, ytemp + 1, itemp]) then
Y_distance[xtemp] := Y_distance[xtemp] + 1
end; {if not last row}
end
else {torus}
begin {note its xtemp, not xtemp+1 now in next line}
if culture^^[xtemp, ytemp, itemp] <> culture^^[(xtemp) mod xmax + 1, ytemp, itemp] then
X_distance[xtemp] := X_distance[xtemp] + 1;
if (culture^^[xtemp, ytemp, itemp] <> culture^^[xtemp, (ytemp) mod ymax + 1, itemp]) then
Y_distance[xtemp] := Y_distance[xtemp] + 1
end;
end; {itemp loop for culture bits}
bit_distance_sum := bit_distance_sum + X_distance[xtemp] + Y_distance[xtemp]; {sum of bit dists at end of periods}
r_dis^^[xtemp, ytemp] := X_distance[xtemp];
d_dis^^[xtemp, ytemp] := Y_distance[xtemp];
if (X_distance[xtemp] = 0) and (xtemp <= Last_X_Boundary) then
n_zero_distance := n_zero_distance + 1;
if (Y_distance[xtemp] = 0) and (ytemp <= Last_Y_Boundary) then
n_zero_distance := n_zero_distance + 1;
end; {xtemp loop for calc distances}
if TextGridReport then
write('y= ', ytemp : 3, '. Dis across: ');
for xtemp := 1 to Last_X_Boundary do
begin {write row for across distances}
if TextGridReport then
write(' ', X_distance[xtemp] : 3);
end; {xtemp for writing row for across distances}
if TextGridReport then
writeln;
if (ytemp <= Last_Y_Boundary) then
if TextGridReport then
write('y= ', ytemp : 3, '. Dis down: ');
if (ytemp <= Last_Y_Boundary) then
begin
for xtemp := 1 to xmax do
begin {write row for distances down}
if TextGridReport then
write(' ', Y_distance[xtemp] : 3);
end; {xtemp for writing row for distances down}
if TextGridReport then
writeln;
end; {if not last row}
end; {ytemp loop}
if TextGridReport then
writeln;
if graphics then
draw_distance_graph;
end;
{ --------------------------------------------------------------- }
function dis (x: x_type; y: y_type; direction: direction_type): integer; {returns bit distance in a direction}
var
xnew, ynew: integer; {to calc wrap-around}
begin {used in Visit}
if flat_map then
begin
if direction = 1 then {dir=North}
dis := d_dis^^[x, y - 1];
if direction = 2 then {dir = E}
dis := r_dis^^[x, y];
if direction = 3 then {dir = W}
dis := r_dis^^[x - 1, y];
if direction = 4 then {dir = S}
dis := d_dis^^[x, y];
end
else {torus. Compare to flat map method above}
begin
if direction = 1 then {dir=North}
begin
if y = 1 then
ynew := ymax
else
ynew := y - 1;
dis := d_dis^^[x, ynew];
end;
if direction = 2 then {dir = E}
dis := r_dis^^[x, y];
if direction = 3 then {dir = W}
begin
if x = 1 then
xnew := xmax
else
xnew := x - 1;
dis := r_dis^^[xnew, y];
{writeln(Region_Output, 'Test333. x = ', x : 2, '. y = ', y : 2, ' ix = ', ix : 2, ' . iy = ', iy : 2, ' . jx = ', jx : 2, ' . jy = ', jy : 2, ' . cult = ', culture^^[jx, jy, 1], culture^^[jx, jy, 2]);}
end;
if direction = 4 then {dir = S}
dis := d_dis^^[x, y];
end;
end;
{ --------------------------------------------------------------- }
procedure Visit (x: x_type; y: y_type; critical_dist: integer); {from Analyze_regions}
var
xq: array[1..10000] of x_type; {space is no bigger than 100x100}
yq: array[1..10000] of y_type;
direction: direction_type;
checkA: boolean;
begin
region_size := 0;
iq := iq + 1; {add node to ready queue}
xq[iq] := x;
yq[iq] := y;
while iq > 0 do {while ready queue not empty...}
begin
ix := xq[iq]; {get node from end of ready queue}
iy := yq[iq];
iq := iq - 1;
for direction := 1 to 4 do {add to ready queue the neighobrs who are legal, not done, and 0 dist }
begin
jx := ix + xmove[direction]; {cacl coords of selected neighbor, j}
jy := iy + ymove[direction];
if not flat_map then {correct for torus}
begin
if jx > xmax then {ok because movement here is only one site}
jx := 1;
if jx = 0 then
jx := xmax;
if jy > ymax then
jy := 1;
if jy = 0 then
jy := ymax;
end;
{write('Test222. ix=', ix : 3, '. iy=', iy : 3, '. jx=', jx : 3, '. jy= ', jy : 3, '. cult=', culture^^[jx, jy, 1], culture^^[jx, jy, 2]);}
{writeln(' done^^[jx,jy ] ', done^^[jx, jy]);}
if (flat_map and (culture^^[jx, jy, 1] <> -1)) then
checkA := true
else
checkA := false;
{writeln(' test 555 ', checkA, flat_map, not flat_map, flat_map = false, checkA or not flat_map, (checkA or (flat_map = false)));}
if ((checkA or not flat_map) and (not done^^[jx, jy])) then
begin
if dis(ix, iy, direction) <= critical_dist then {crit distance = 0 for region, bitmax-1 for zone}
begin
iq := iq + 1; {add this valid neighbor to ready queue}
xq[iq] := jx;
yq[iq] := jy;
done^^[jx, jy] := true; {change its status to done}
region_size := region_size + 1;
{writeln('test 11/21. jx, jy, culture^^[jx, jy, 1]', jx, jy, culture^^[jx, jy, 1]);}
end; {if dis< min}
end;{if legal and not done}
end;{direction}
end; {while iq>0}
if region_size = 0 then
region_size := 1; {needed for 1x1 regions which have no valid neighbors}
end;
{ --------------------------------------------------------------- }
procedure Ck_Equal_Cultures; {from Analyze Regions}
var
x, x2: x_type; {}
y, y2: y_type;
bit: integer; {local index of cultural bit}
r, r2: integer; {local indices of regions}
cultures_same: boolean; {two regions have same culture}
begin
for r := 1 to regions do
begin
x := x_region^^[r];
y := y_region^^[r];
for r2 := 1 to r - 1 do
begin
x2 := x_region^^[r2];
y2 := y_region^^[r2];
bit := 0;
cultures_same := true;
repeat
begin
bit := bit + 1;
if culture^^[x, y, bit] <> culture^^[x2, y2, bit] then
cultures_same := false;
end;
until (not cultures_same) or (bit = bitmax);
if cultures_same then
{writeln('Culture of (', x : 2, y : 2, ')=Culture of (', x2 : 2, y2 : 2, ').'); }
end;
end;
end;
{ --------------------------------------------------------------- }
procedure Report_Per_Region (x: x_type; y: y_type); {from Analyze_Regions}
var
r2: integer; {regional count}
x2: x_type;
y2: y_type;
culture_dis: integer;
bit_temp: bit_type;
begin
write('R ', regions : 3, ' (', x : 2, ',', y : 2, '). Sz =', region_size : 3);
write('. Cul= ');
for bit_temp := 1 to bitmax do
begin
write(culture^^[x, y, bit_temp] : 2);
end;
write('. D to prev regs = ');
for r2 := 1 to regions - 1 do
begin
x2 := x_region^^[r2];
y2 := y_region^^[r2];
culture_dis := 0;
for bit_temp := 1 to bitmax do
begin
if culture^^[x, y, bit_temp] <> culture^^[x2, y2, bit_temp] then
culture_dis := culture_dis + 1;
end;
write(culture_dis : 3);
end;{r}
writeln;
end;
{ --------------------------------------------------------------- }
function Map_Clear: boolean; {are all boundaries 0 or bitmax?}
{not used starting Org 14}
var
boundaries, distance_if_clear: longint;
begin
boundaries := (Last_X_Boundary) * ymax + xmax * (Last_Y_Boundary); {# boundaries down and across}
distance_if_clear := bitmax * (boundaries - n_zero_distance);
if distance_if_clear = bit_distance_sum then
map_clear := true
else
map_clear := false;
end;
{ --------------------------------------------------------------- }
function Neighbor_Clear: boolean; {are all neighobrs 0 or bitmax?}
var {added at Org14}
ix, iy, jx, jy, direction, bit: integer; {temp local variables}
match_count: integer;
label
1;
begin
neighbor_clear := true; {assume all neighbors are 0 or bitmax from all sites}
for ix := 1 to Xmax do
begin
for iy := 1 to Ymax do
begin
for direction := 1 to Interaction_Types do
begin
if flat_map then
begin
jx := ix + xmove[direction]; {cacl coords of selected neighbor, j}
jy := iy + ymove[direction]
end
else {torus map}
begin {note: have to add Xmax or Ymax to be sure to get positive number}
jx := (ix + xmove[direction] - 1 + Xmax) mod xmax + 1; {cacl coords of selected neighbor, j}
jy := (iy + ymove[direction] - 1 + Ymax) mod ymax + 1;
end;
if ((flat_map = false) or ((jx > 0) and (jx < XMaxPlus1) and (jy > 0) and (jy < YMaxPlus1))) then
begin {legal neighbor. Note any neigh is legal if flat_map= false }
match_count := 0;
for bit := 1 to bitmax do
begin
if (culture^^[ix, iy, bit] = culture^^[jx, jy, bit]) then
match_count := match_count + 1;
end; {bit}
if (match_count > 0) and (match_count < bitmax) then
begin
neighbor_clear := false;
goto 1; {leave if a neighbor is not clear}
end;
end; {legal neighbor}
end;{neighbor}
end; {it}
end; {ix}
1: {go here here if a neighbor is not clear}
end;
{ --------------------------------------------------------------- }
procedure Report_Final_Region (x: x_type; y: y_type); {from Analyze_Regions}
var {done if (WriteRegionStatsAtEnd) and (neighbor_clear = true)}
r2: integer; {regional count}
x2: x_type;
y2: y_type;
culture_dis: integer;
bit_temp: bit_type;
begin
write(region_output, case_number : 3, ' ', map : 3, ' ', pop : 3);
write(region_output, ' ', regions : 3, ' ', x : 2, ' ', y : 2, ' ', region_size : 3);
for bit_temp := 1 to bitmax do
begin
write(region_output, ' ', culture^^[x, y, bit_temp] : 2);
end;
for r2 := 1 to regions - 1 do
begin
x2 := x_region^^[r2];
y2 := y_region^^[r2];
culture_dis := 0;
for bit_temp := 1 to bitmax do
begin
if culture^^[x, y, bit_temp] <> culture^^[x2, y2, bit_temp] then
culture_dis := culture_dis + 1;
end;
write(region_output, ' ', culture_dis : 3);
end;{r}
writeln(region_output);
end;
{ --------------------------------------------------------------- }
procedure Analyze_regions; {count number of regions, and types}
var {based on breadth first search.See Stubbs and Webre, Data Structures, 359ff}
x: x_type;
y: y_type;
begin
regions := 0;
iq := 0; {initialize ready queue}
for x := 1 to xmax do
begin
for y := 1 to ymax do
begin
done^^[x, y] := false;
end;
end; {end initialization}
for x := 1 to xmax do {for each node...}
begin
for y := 1 to ymax do
begin
if not done^^[x, y] then {if not done then record it and visit it}
begin
regions := regions + 1; {count regions}
Visit(x, y, 0); {0 is min distance allowed for regions}
if WriteRegionStats then
Report_Per_Region(x, y);
if (WriteRegionStatsAtEnd) and (neighbor_clear = true) then
Report_Final_Region(x, y);
x_region^^[regions] := x; {add loc to list of regions}
y_region^^[regions] := y;
end;
end;
end; {for each node}
Ck_Equal_Cultures;
{writeln('Period ', period : 4, '. Number of regions = ', regions : 4);}
end;
{ --------------------------------------------------------------- }
procedure Analyze_zones; {count number of zones, and types. mod. from Analyze regions}
var {based on breadth first search.See Stubbs and Webre, Data Structures, 359ff}
x: x_type;
y: y_type;
max_distance: integer;
begin
zones := 0;
iq := 0; {initialize ready queue}
max_distance := bitmax - 1; {a zone may have any distance, except maxmal}
for x := 1 to xmax do
begin
for y := 1 to ymax do
begin
done^^[x, y] := false;
end;
end; {end initialization}
for x := 1 to xmax do {for each node...}
begin
for y := 1 to ymax do
begin
if not done^^[x, y] then {if not done then record it and visit it}
begin
zones := zones + 1; {count regions}
Visit(x, y, max_distance); {0 is min distance allowed for regions}
x_region^^[regions] := x; {add loc to list of regions}
y_region^^[regions] := y;
end;
end;
end; {for each node}
Ck_Equal_Cultures;
{writeln('Period ', period : 4, '. Number of regions = ', regions : 4);}
end;
{ --------------------------------------------------------------- }
procedure Output_period; {output for a single period}
var
total_cycles: longint;
begin
if TextGridReport then
report_output_period_test; {test output for this perod}
CalcDistance; {Calc distance between neighbors}
Analyze_regions; {Analyze number of regions, and their cultures}
if WritePeriodicStats then
begin
Analyze_zones; {Analyze number of zones - ie totally incompatible domains}
total_cycles := period * cyclemax;
Write(period_output, case_number : 5, ' ', map : 5, ' ', pop : 5, ' ', period : 5, ' ', total_cycles : 8);
Write(period_output, ' ', bit_distance_sum : 5, ' ', n_zero_distance : 5, ' ', regions : 3, ' ', changes_this_period : 5);
Writeln(period_output, ' ', zones : 3);
end;
if button then
begin
writeln('**** Halted by user ****');
writeln(final_output, '**** Halted by user ****');
halt;
end;
end;
{ --------------------------------------------------------------- }
procedure Output_pop; {output for a single population}
var
x, y: integer;
begin
{write('End of Pop ', pop : 3, '. Per ', period : 3, ' . Map_Clear =', Map_Clear : 2, '. ', regions : 3, ' regions. ');}
{writeln(bit_distance_sum : 4, ' total dis. ', n_zero_distance : 4, ' # zero dis.');}
write(Case_Number : 4);
if mapmax > 1 then
write(' ', map : 3);
write(' ', pop : 5, ' ', Xmax : 4, ' ', Interaction_types : 4, ' ', bitmax : 4, ' ', allelemax : 4, ' ');
write(Map_Clear : 5, ' ', period : 6, ' ', Xmax : 4, ' ', regions, ' ');
writeln(bit_distance_sum : 8, ' ', n_zero_distance : 4);
write(final_output, Case_Number : 4);
if mapmax > 1 then
write(final_output, ' ', map : 3);
write(final_output, ' ', pop : 5, ' ', Xmax : 4, ' ', Interaction_types : 4, ' ', bitmax : 4, ' ', allelemax : 4, ' ');
write(final_output, Map_Clear : 5, ' ', period : 6, ' ', Xmax : 4, ' ', regions, ' ');
writeln(final_output, bit_distance_sum : 8, ' ', n_zero_distance : 4);
if report_bit_origin then {report first bit's origin}
begin
for y := 1 to ymax do
begin
write('Origin of first bit, y= ', y : 3, '. ');
for x := 1 to xmax do
begin
write(origin_first_bit^^[x, y] : 5);
end; {x}
writeln;
end;{y}
end; {report_bit_origin}
{ Write('TEST: String Length = ', bitmax : 3, '. Alleles=', allelemax : 3); }
{ Write('. Prob diff=(1-1/A)**L=', exp(bitmax * ln(1 - 1 / allelemax)) : 10 : 4);}
{ Writeln('. Reciprocal =', 1 / exp(bitmax * ln(1 - 1 / allelemax)) : 10 : 4); }
{ writeln(' Test . Pop ', pop : 3, ' . Culture [ 011 ] = ', culture^^[0, 1, 1] : 2, ' . Culture [ 341 ] = ', culture^^[3, 4, 1] : 2);}
end;
{ --------------------------------------------------------------- }
procedure initialize_edge_of_map; {enter -1 just beyond edges; for regional analysis}
var
x, y, bit: integer;
begin
for x := 1 to xmax do {initialize ind's beyond the internal borders. Never changes.}
begin
for bit := 1 to bitmax do
begin
culture^^[x, 0, bit] := -1; {initial culture for beyond top internal border}
culture^^[x, ymax + 1, bit] := -1; {initial culture for beyond bottom internal }
end;
end; {x}
for y := 1 to ymax do {initialize indiv's beyond the internal borders. Never changes.}
begin
for bit := 1 to bitmax do
begin
culture^^[0, y, bit] := -1; {initial culture for beyond left internal border}
culture^^[Xmax + 1, y, bit] := -1; {initial for beyond right internal border}
end;
end;{y}
end;
{ --------------------------------------------------------------- }
procedure Output_Run; {Output for run - last thing written}
begin
gettime(end_datetime);
end_hour := end_datetime.hour; {to force long interger}
end_time := 60 * 60 * end_hour + 60 * end_datetime.minute + end_datetime.second;
duration := end_time - start_time;
Writeln('Duration of this run is ', duration : 5, ' seconds.');
Writeln(final_output, 'Duration of this run is ', duration : 5, ' seconds.');
{Writeln('Broadcast count=', broadcast_count : 4); }
end;
{ --------------------------------------------------------------- }
procedure Cycle_Loop; {Assume Mthod B: do only one interaction when there is match}
var { Allow broadcasting}
cycle: longint; {current cycle (ie active indiv), 1..cyclemax}
begin
for cycle := 1 to cyclemax do {CYCLE of one active actor}
begin
ix := Random_one_to_n(Xmax); {select X coord of active actor}
iy := Random_one_to_n(Ymax); {select Y coord of active actor}
bit := random_one_to_n(bitmax); {bit to be checked for match}
if broadcast_allowed then {do broadcast interaction}
if random < broadcast_cutoff_number then
interact_from_broadcast
else
interact_from_neighbor {braodcast allowed, but didn't happen this time}
else {do interaction with neighbor}
interact_from_neighbor; {broadcast not allowed}
if mutuation_rate_perpop_percycle > 0.0 then
mutate; {mutate culture bits for drift}
end; { cycle of one active }
end;
{- - - - - - - - - - - - - - - - - ----------------- - - - - - - - - - - - - - - - - }
{ --------------------------------------------------------------- }
{M A I N P R O G R A M }
begin
initialize_run; {initialize whole run}
case_number := 0; {sequential count of pops}
for Xmax_index := 1 to N_Xmax do
begin
Xmax := Xmax_value[Xmax_index];
Ymax := Xmax; {assume square map}
XmaxPlus1 := Xmax + 1;
YmaxPlus1 := Ymax + 1;
XmaxYmax := Xmax * Ymax;
if flat_map then
begin
Last_X_Boundary := Xmax - 1; {flat map}
Last_Y_Boundary := Ymax - 1;
end
else
begin
Last_X_Boundary := Xmax; {torus}
Last_Y_Boundary := Ymax;
end;
if graphics then
frame_rectangle;
for interaction_types_index := 1 to N_interaction_types do
begin
interaction_types := interaction_types_value[interaction_types_index];
for bitmax_index := 1 to N_bitmax do
begin
bitmax := bitmax_value[bitmax_index];
initialize_edge_of_map;
for allelemax_index := 1 to N_allelemax do
begin
allelemax := allelemax_value[allelemax_index];
for map := 1 to mapmax do {generate map starts}
map_seed[map] := random;
for pop := 1 to popmax do {generate interaction starts}
interaction_seed[pop] := random;
for map := 1 to mapmax do {do a map}
begin
for pop := 1 to popmax do {do a pop with current map}
begin
randseed := map_seed[map];
{writeln('Map Repeat Test: randseed for this map=', randseed : 8);}
Initialize_pop; {INITIALIZE current population}
randseed := interaction_seed[pop]; {to start interactions}
{writeln('Interaction Repeat Test: randseed to start interactions=', randseed : 8);}
case_number := case_number + 1; {sequential count of pops's}
period := 0;
repeat {REPORTING PERIOD }
begin
period := period + 1;
Initialize_Period;
Cycle_Loop;
Output_Period; {output of a period, incl regional analysis}
end;{period}
{until (period = periodmax) or (neighbor_clear = true); }
until (period = periodmax); {changed from prev line 10/30/96 for mutation}
{end pop at periodmax or neighbors are all 0 or bitmax}
Output_Pop; {output of the population}
end; {pop}
end; {map}
end; {bitmax_index}
end;{allelemax_index}
end;{interaction_type_index}
end; {Xmax_index}
Output_Run; {output of the run}
end.{main program}
University of Michigan Center for the Study of Complex Systems
Contact cscs@umich.edu.
Revised December 20, 1996.