* Initially there is X number of people on a map and a few werewolves.
* Werewolves head toward whatever unprotected humans are on the map.
* If they encounter a human they kill that human and go towards they next unprotected human.
* Humans do not know about the werewolves until they become within sensory range - specified by the user.
* Once they become within sensory range humans immediately seek the closest shelter.
* Shelters are distributed throughout the map.
* Once a human is in a shelter the human is safe.
***************************************************************
* Simulation Parameters Begin
***************************************************************
clear
* Set the number of humans
local num_humans=1500* Human sight range (radius)
local human_vision=65
* Human run speed
local human_speed = 1
* Set number of shelters
local num_shelters = 3
* Specify the initial number of werewolves
local num_werewolves = 8
* Werewolf run speed (if humans are faster than werewolf then it is only the unlucky human that will ever be caught)
local werewolf_speed = 5
* Specify the dimensions of the world grid
local xrange=400
local yrange=400
* Graph interval every x rounds the program graphs what happens
local graph_int = 20
* Tell stata to not graph as the simulation progresses
local draw=1
***************************************************************
* Simulation Parameters End
***************************************************************
* This tells stata to either draw the graphs as it runs the code or wait till the end to draw the graphs.
if `draw'==0 local nodraw nodraw
set obs `num_humans'
gen hmn_id=_n
label var hmn_id "Human ID"
* To begin with all types are human
gen ty = 1
lab var ty "Observation type 1=human, 2=werewolf, 3=shelter"
label define obj_type 1 "Human" 2 "Werewolf" 3 "Shelter" 4 ///
"Deceased" 5 "Spooked Humans" 6 "Safe"
label values ty obj_type
* Then we will add the werewolves
set obs `=`num_humans'+`num_werewolves''
replace ty = 2 if ty==.
* Put the werewolves on the top of the list
gsort -ty
gen were_id = _n if ty==2
label var were_id "Werewolf ID"
* Next the shelters
set obs `=`num_humans'+`num_werewolves'+`num_shelters''
replace ty = 3 if ty==.
tab ty
gsort -ty
gen shelter_id = _n if ty==3
label var shelter_id "Shelter ID"
* Looks pretty good so far.
* People shelters and werewolves are uniformly throughout the data.
gen x=runiform()*`xrange'
gen y=runiform()*`yrange'
* Though werewolves will enter the maps from the edges
gen xy_edge = rbinomial(1,.5)
replace x=round(x, `xrange') if xy_edge==0 & ty == 2
replace y=round(y, `yrange') if xy_edge==1 & ty == 2
* This makes it so the werewolves start on one of the sides of the map.
* Generate
sum x
replace x=rnormal(r(mean),r(sd)*5/8) if ty == 3
sum y
replace y=rnormal(r(mean),r(sd)*5/8) if ty == 3
two (scatter y x if ty==1 , mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==3, msymbol(square) color(yellow)) ///
(scatter y x if ty==2, color(cranberry) ) ///
, legend(label (1 "Humans") ///
label (2 "Werewolves") ///
label (3 "Shelters") ///
rows(1)) `nodraw' title(Full Moon - Round 0) ///
plotregion(fcolor(gs2))
***************************************************************
* First we will calculate which direction and where the humans will run.
* This direction is constant since humans do not try to avoid werewolves, they just run for shelter.
***************************************************************
* We will loops through shelters and calculate distance for each human.
global shelter_distance_list
qui forv i=1(1)`num_shelters' {
* The first x number of observations are the werewolves.
sort shelter_id
local shelter_x = x[`i']
local shelter_y = y[`i']
* were_x and were_y hold the position of the current werewolf looking for its prey
* Now let's calculate the distance of the spooked humans from the shelters.
cap drop dist_hm_shltr`i'
gen dist_hm_shltr`i'=((x-`shelter_x')^2 + (y-`shelter_y')^2)^.5 if ty==1
global shelter_distance_list $shelter_distance_list , dist_hm_shltr`i'
}
di "$shelter_distance_list"
* We can see in out list of shelters there is unwanted comma.
* For now we will just fill in an meaningless value to prevent error.
* Now we figure out which shelter is the closest for all of the humans collectively
gen closest_shelter = .
gen closest_shelter_x = .
gen closest_shelter_y = .
gen closest_shelter_dist = .
qui forv i=1(1)`num_shelters' {
replace closest_shelter = `i' if dist_hm_shltr`i' == min(99999 $shelter_distance_list)
replace closest_shelter_dist = dist_hm_shltr`i' if closest_shelter == `i'
local shelter_x = x[`i']
local shelter_y = y[`i']
replace closest_shelter_x = `shelter_x' if closest_shelter == `i'
replace closest_shelter_y = `shelter_y' if closest_shelter == `i'
}
tab closest_shelter
* Now we will calculate the ideal trajectory for each human given that human is running for shelter.
* We know a, b, c [see graph] and the next section for algebra.
gen a = y-closest_shelter_y
gen b = x-closest_shelter_x
gen a_prime=1/(1+b^2/a^2)^.5
* b'^2 = 1-a'^2
* b' = (1-a'^2)^.5
gen b_prime = (1-(a_prime)^2)^.5
gen x2 = x - b_prime*sign(b) if ty == 1
gen y2 = y - a_prime*sign(a) if ty == 1
replace x2 = x2 - b_prime*sign(b) if ty == 1
replace y2 = y2 - a_prime*sign(a) if ty == 1
two (scatter y x if ty==1 , mcolor(gs12) msymbol(smcircle)) ///
(scatter y2 x2 if ty==1, mcolor(sandb) msymbol(smcircle)) ///
(scatter y x if ty==3, msymbol(square) color(yellow)) ///
, legend(label (1 "Starting Place") ///
label (2 "One step") ///
label (3 "Shelters") ///
rows(1)) `nodraw' title(Sample humans running for shelter) ///
legend(color(white) region(fcolor(black))) ///
plotregion(fcolor(gs2))
drop x2 y2
***************************************************************
* Werewolves search for human targets
***************************************************************
* I loop through each of the werewolves.
gen x2=.
gen y2=.
* Each werewolf will survey all of the humans and find the one who is the closest.
qui forv i=1(1)`num_werewolves' {
* The first x number of observations are the werewolves.
sort were_id
local were_x = x[`i']
local were_y = y[`i']
* were_x and were_y hold the position of the current werewolf looking for its prey
* Now let's calculate the distance of that werewolf from the humans.
cap drop dist_hm_wr
gen dist_hm_wr`i'=((x-`were_x')^2 + (y-`were_y')^2)^.5 if ty==1 | ty==5
* We take the square but for finding targets it is unnecessary since the ranks of the different distances remains constant.
replace ty=5 if dist_hm_wr`i' < `human_vision'
sort dist_hm_wr
* The top of this list is the target human.
* If a human is less than `werewolf_speed' distance from the were and it is the closest to that were then the human is dead.
replace ty = 4 if dist_hm_wr < `werewolf_speed' & _n==1
* Now we need to calculate the attack vector.
* Ie, what is the linear projection that gets the werewolf from its current location to the human fastest (in terms of x and y movement)?
* We can use some algebra and a little trig to find the proper vector.
* We know a, b, c
local a = `were_y'-y[1]
local b = `were_x'-x[1]
* And c' which is equal to the speed of werewolves
local c_prime=`werewolf_speed'
* Triangle abc is similar to triangle a'b'c' so the ratios b/a=b'/a'
* b'=b(a'/a)
* (1) b'^2=b^2 * a'^2 / a^2
* (2) Pathagoras (500BC): a'^2 + b'^2 = c'^2
* (3) instert (1) into (2): a'^2 + b^2 * a'^2 / a^2 = c'^2
* a'^2(1 + b^2 / a^2) = c'^2
* a'^2 = c'^2/(1 + b^2 / a^2)* a'=c'/(1 + b^2 / a^2)^.5
local a_prime=`c_prime'/(1 + (`b')^2 / (`a')^2)^.5
* b'^2 = c'^2-a'^2
* b' = (c'^2-a'^2)^.5
local b_prime = ((`c_prime')^2-(`a_prime')^2)^.5
* It is important to note the direction that the werewolf needs to travel as well.
/*
replace x = x + `b_prime'*sign(`b') if were_id == `i'
replace y = y + `a_prime'*sign(`a') if were_id == `i'
*/
replace x2 = x - `b_prime'*sign(`b') if were_id == `i'
replace y2 = y - `a_prime'*sign(`a') if were_id == `i'
noi di "Werewolf `i' moves"
}
two (scatter y x if ty==1 , mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==5 , color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry)) ///
(scatter y2 x2 if ty==2, color(red)) ///
(scatter y x if ty==3, msymbol(square) color(yellow)) ///
(scatter y x if ty==4, msymbol(smplus) color(gs16)) ///
, legend(label (1 "Unconcerned Humans") ///
label (2 "Spooked Humans") ///
label (3 "Werewolves") ///
label (4 "Werewolf's 1st Move") ///
label (5 "Shelters") ///
label (6 "Deceased") ///
rows(2)) `nodraw' title(Full Moon - Round 0) ///
legend(color(white) region(fcolor(black))) ///
plotregion(fcolor(gs2))
**************************************************************** Spooked Humans Seek Out Closest Shelter
***************************************************************
* Move spooked humans towards shelters.
replace closest_shelter_dist = closest_shelter_dist-`human_speed'
replace ty=6 if closest_shelter_dist < 0
* Once a human gets to a shelter, that human is considered safe.
replace x2 = x - b_prime*sign(b)*`human_speed' if ty == 5
replace y2 = y - a_prime*sign(a)*`human_speed' if ty == 5
two (scatter y x if ty==1 | ty==5, mcolor(gs12) msymbol(smcircle)) ///
(scatter y2 x2 if ty==5 , color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry)) ///
(scatter y2 x2 if ty==2, color(red)) ///
(scatter y x if ty==3, msymbol(square) color(yellow)) ///
(scatter y x if ty==4, msymbol(smplus) color(gs16)) ///
, legend(label (1 "Unconcerned Humans") ///
label (2 "Spooked Humans") ///
label (3 "Werewolves") ///
label (4 "Werewolf's 1st Move") ///
label (5 "Shelters") ///
label (6 "Deceased") ///
rows(2)) `nodraw' title(Full Moon - Round 0) ///
legend(color(white) region(fcolor(black))) ///
plotregion(fcolor(gs2))
**** LOOKS like things are working properly
* Let's make the position changes to x and y permanent
replace x = x2 if ty == 5 | ty == 2
replace y = y2 if ty == 5 | ty == 2
drop x2 y2
***************************************************************
* Let's allow for some summary statistics
gen obs_num = _n
gen unaware_humans = `num_humans' if _n == 1
gen spooked_humans = 0 if _n == 1
gen deceased_humans = 0 if _n == 1
gen safe_human = 0 if _n == 1
***************************************************************
* Now we will begin the system loop
***************************************************************
global graph_list
cap
local ii=0
while _rc==0 {
local ii=`ii'+1
noi di "Round `=`ii'+1'"
***************************************************************
* Let's record some summary statistics
qui sum ty if ty==1
replace unaware_humans = r(N) if obs_num == `ii'+1
qui sum ty if ty==5
replace spooked_humans = r(N) if obs_num == `ii'+1
qui sum ty if ty==4
replace deceased_humans = r(N) if obs_num == `ii'+1
qui sum ty if ty==6
replace safe_human = r(N) if obs_num == `ii'+1
***************************************************************
* First the werewolf moves
* Each werewolf will survey all of the humans and find the one who is the closest.
sum ty if ty==1 | ty==5
if r(N)==0 cap end_loop
if r(N)>0 qui forv i=1(1)`num_werewolves' {
* The first x number of observations are the werewolves.
sort were_id
local were_x = x[`i']
local were_y = y[`i']
* were_x and were_y hold the position of the current werewolf looking for its prey
* Now let's calculate the distance of that werewolf from the humans.
cap drop dist_hm_wr
gen dist_hm_wr`i'=((x-`were_x')^2 + (y-`were_y')^2)^.5 if ty==1 | ty==5
* We take the square but for finding targets it is unnecessary since the ranks of the different distances remains constant.
replace ty=5 if dist_hm_wr`i' < `human_vision'
sort dist_hm_wr
* The top of this list is the target human.
* If a human is less than `werewolf_speed' distance from the were and it is the closest to that were then the human is dead.
replace ty = 4 if dist_hm_wr < `werewolf_speed' & _n==1
* Now we need to calculate the attack vector.
* Ie, what is the linear projection that gets the werewolf from its current location to the human fastest (in terms of x and y movement)?
* We can use some algebra and a little trig to find the proper vector.
* We know a, b, c
local a = `were_y'-y[1]
local b = `were_x'-x[1]
* And c' which is equal to the speed of werewolves.
local c_prime=`werewolf_speed'
* Triangle abc is similar to triangle a'b'c' so the ratios b/a=b'/a'
* b'=b(a'/a)
* (1) b'^2=b^2 * a'^2 / a^2
* (2) Pathagoras (500BC): a'^2 + b'^2 = c'^2
* (3) insert (1) into (2): a'^2 + b^2 * a'^2 / a^2 = c'^2
* a'^2(1 + b^2 / a^2) = c'^2
* a'^2 = c'^2/(1 + b^2 / a^2)
* a'=c'/(1 + b^2 / a^2)^.5
local a_prime=`c_prime'/(1 + (`b')^2 / (`a')^2)^.5
* b'^2 = c'^2-a'^2
* b' = (c'^2-a'^2)^.5
local b_prime = ((`c_prime')^2-(`a_prime')^2)^.5
replace x = x - `b_prime'*sign(`b') if were_id == `i'
replace y = y - `a_prime'*sign(`a') if were_id == `i'
}
***************************************************************
* Then the humans move
* Move spooked humans towards shelters.
replace closest_shelter_dist = closest_shelter_dist-`human_speed' if ty == 5
replace ty=6 if closest_shelter_dist < 0 & ty == 5
* Once a human gets to a shelter, that human is considered safe.
replace x = x - b_prime*sign(b)*`human_speed' if ty == 5
replace y = y - a_prime*sign(a)*`human_speed' if ty == 5
***************************************************************
* Then the humans move
if `ii' / `graph_int'==int(`ii' / `graph_int') {
two (scatter y x if ty==1, mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==5, color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry)) ///
(scatter y x if ty==3, msymbol(square) color(yellow)) ///
(scatter y x if ty==4, msymbol(smplus) color(gs16)) ///
, legend(label (1 "Unconcerned Humans") ///
label (2 "Spooked Humans") ///
label (3 "Werewolves") ///
label (4 "Shelters") ///
label (5 "Deceased") ///
rows(2)) `nodraw' title(Full Moon - Round `ii') ///
legend(color(white) region(fcolor(black))) ///
plotregion(fcolor(gs2))
two (scatter y x if ty==1, mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==5, color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry)) ///
(scatter y x if ty==3, msymbol(square) color(yellow)) ///
(scatter y x if ty==4, msymbol(smplus) color(gs16)) ///
, legend(off) nodraw name(round`ii', replace) ///
plotregion(fcolor(gs2))
global graph_list $graph_list round`ii'
}
}
***************************************************************
* Now we will end the system loop
***************************************************************
sort obs_num
two (line unaware_humans round, sort) ///
(line spooked_humans round, sort) ///
(line deceased_humans round, sort) ///
(line safe_human round, sort) ///
, legend(label (1 "Unconcerned Humans") ///
label (2 "Spooked Humans") ///
label (3 "Deceased Humans") ///
label (4 "Safe Humans"))
* Different specifications at the beginning of the model yield different results. (werewolves always win though)
sleep 5000
graph combine $graph_list* (c) Francis Smart 2012.
No comments:
Post a Comment