Thursday, June 14, 2012

Werewolf attack: Spatial Multi-agent simulation - basic artificial intelligence

* 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.

No comments:

Post a Comment