* This is a simulation of a werewolf attack using Stata.

* 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

gen round=obs_num if unaware_humans!=.

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.