A Sudoku Puzzle Solver - attempt 1

I have programmed up a R based Sudoku problem solver for Sudoku puzzles of that only require simple inference.  In these puzzles a solution can be found using only first order inference.  This solver can be found at the end of the code located in the git:

https://github.com/EconometricsBySimulation/2013-06-14-Sudoku

A puzzle of this form can look like this (generated from the R code on the git)

A solution to this problem can be found using the following steps generated from the solver:

[1] "Round 1 -sub out 1 4 with 4"  That is sub out row 1 column 4 with a 4.

[1] "Round 1 -sub out 1 6 with 5"

[1] "Round 1 -sub out 2 7 with 5"

[1] "Round 1 -sub out 2 8 with 3"

[1] "Round 1 -sub out 2 9 with 8"

[1] "Round 1 -sub out 3 5 with 9"

[1] "Round 1 -sub out 3 9 with 2"

[1] "Round 1 -sub out 4 7 with 2"

[1] "Round 1 -sub out 5 5 with 4"

[1] "Round 1 -sub out 5 8 with 6"

[1] "Round 1 -sub out 6 1 with 2"

[1] "Round 1 -sub out 6 6 with 9"

[1] "Round 1 -sub out 7 5 with 5"

[1] "Round 1 -sub out 8 9 with 9"

[1] "Round 1 -sub out 9 1 with 3"

[1] "Round 1 -sub out 9 2 with 2"

[1] "Round 1 -sub out 9 3 with 4"

[1] "Round 1 -sub out 9 6 with 1"

[1] "Round 1 -sub out 9 9 with 7"

[1] "Round 2 -sub out 1 1 with 1"

[1] "Round 2 -sub out 1 2 with 3"

[1] "Round 2 -sub out 1 3 with 2"

[1] "Round 2 -sub out 2 1 with 7"

[1] "Round 2 -sub out 3 1 with 6"

[1] "Round 2 -sub out 3 2 with 5"

[1] "Round 2 -sub out 3 8 with 1"

[1] "Round 2 -sub out 4 2 with 8"

[1] "Round 2 -sub out 4 4 with 6"

[1] "Round 2 -sub out 4 8 with 9"

[1] "Round 2 -sub out 4 9 with 5"

[1] "Round 2 -sub out 5 2 with 7"

[1] "Round 2 -sub out 5 3 with 5"

[1] "Round 2 -sub out 5 4 with 2"

[1] "Round 2 -sub out 5 9 with 3"

[1] "Round 2 -sub out 7 1 with 8"

[1] "Round 2 -sub out 7 2 with 9"

[1] "Round 2 -sub out 7 3 with 6"

[1] "Round 2 -sub out 8 4 with 8"

The results look like this.  You can see that it was not able to solve all of the missing blocks though it is pretty close.

Interestingly, I think this puzzle has two potential solutions.

This solver is not at the level I wanted it to be at.  I can imagine two more levels of inference that it could handle: one in which the solution uses secondary inference information.  Such as knowing that if it chooses a particular value the other block cannot be filled with that value.  This I call a type 2 and finally a guess and check method which I call a type 3.  The guess and check method should be easy to program but have the difficulty that if there are many squares not solved by the previous methods, it will potentially bog down the machine as it has to check a near infinite number of potential solutions.

However, there should be very few puzzles in which guess and check would be required since those puzzles would be even harder for most people to solve than they would for an algorithm.

I am not satisfied with this solution but I wanted to post this because I had worked it out a few weeks ago and my real work is calling so I am not sure when I will get back to this.

Sudoku Automation Solver Challenge - R

On a recent flight I was bored waiting for the plane to land and I tried out the electronic sudoku game that they had offered.  I found the game surprisingly interesting as I realized that it is far more entertaining when you cannot use paper or pencil to augment the cognitive solution seeking process.  In addition, having your score based on your time presented an objective challenge and measure to test yourself against.

Figure 1:a sample Sudoku grid

After that I decided to design a Sudoku generator.  Which I have included:

https://github.com/EconometricsBySimulation/2013-06-14-Sudoku

It has been an entertaining exercise figuring out methods to produce matrices that conform to the specifications of a Sudoku table.  That is a table which distributes sets of 1 to 9 distributed both horizontally, vertically, and in a 3x3 grid without duplicates for all squares in a 9x9 grid.

Figure 2: Ultimate Sudoku

You can see from the R script (and the grid produced by it figure 1) that I am able to accomplish this very easily through a randomized guess and check process (with the computer finding a solution in 1 in 50-500 attempts) that is capable of drawing nearly an infinite number of grids (of course there is a finite though very large number of possible grids).

I have also made a variant that I call Ultimate Sudoku which conforms to Sudoku rules and also adds one more rule to the grid in which each position in a 3x3 set is part of a set of all other 3x3 sets in which there can be no duplicate entries from the set.  That is for instance that the first square in a 3x3 grid cannot have duplicate entries with any of the other first squares in the other grid.  An example table that I generated using R is figure 2.

I don't think it is likely that Ultimate Sudoku will take off but for those hard core Sudoku players this variate would present an option that would allow for more challenging games since players must keep track of four overlapping sets of elements rather than 3.  Under this set up effectively more spaces could be made blank.  The R script is also capable of generating these though it is a bit more challenging to find grids that conform to all four rules (usually between 1 out of 4-10 thousand attempts).

I have also defined a function to strategically remove values from the plotted grid so that players are faced with an incomplete grid.  I have been able to vary the difficulty of these grids as well by specifying the average number of missing elements (from the set of 1 through 9) for each grid value.  You can see examples of problem grids I have generated in figures 3 and 4.

Figures 3 and 4

What I realized though was what my system is still not perfect because it is possible that I can come up with grids which might not have solutions which would present a problem.  What is really needed is an algorithm which is capable of doing what Sudoku players are capable of doing, that is solve the tables.  So rather that taking the next step designing an algorithm to solve my sample tables I thought that perhaps some of my readers out there might be interested in attempting to design an algorithm to solve these potential Sudoku tables and identify which ones are not feasible.

So, that is the challenge!  It should be easy to run the code an generate and map the tables using the Git.

Spatial Critter Swarming Simulation

# I am interested in how small bits of individualized instructions can create collective action.

# In this simulation I will give a single instruction to each individual in the swarm.

# Choose another individual who is not too close, then accelerate towards that individual.

# I also control momentum causing the previous movement and direction to only decay at a small rate.

# TO SEE Original Script


# Critters are initially distributed randomly on a 1 x 1 grid.

ncritters = 40

xypos = matrix(runif(ncritters*2),ncol=2)
plot(xypos, main="Critters are Initially Distributed Randomly"
          , xlab="X", ylab="Y")

# Now let's imagine that each critter has an ideal safe distance from each other critter.

safe.dist = .3

critter.speed = .001

# If another critter is not at that safe distance than the critter will move towards the closest nearby critter.

# Let's see how this works.

# First let's check how close each critter is to each other critter.
# We will accomplish this by going through each critter and checking how far away each other critter is.
distances = NULL
for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

# In order to prevent critters from always chasing whatever is closest to them (and themselves) we drop anything which is closer than the safe.distance.
distances[abs(distances)closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]
  # The apply command will apply the order command to each row whiel the [1,] selects only the critter that is closes.

# Plot the
plot(xypos, xlab = "X", ylab = "Y")
for (i in 1:ncritters) arrows(x0=xypos[i,1], y0=xypos[i,2],
                              x1=xypos[closest,][i,1],
                              y1=xypos[closest,][i,2],
                              length=.1)

# This calculates the difference between the current position of each critter and that of the closest critter.
ab = xypos-xypos[closest,]

# To see how this is calculated, see my previous post simulating a werewolf attack.

# Now calculate the difference in the horizontal and vertical axes that the critters will move as a projection into the direction of the closest critter outside of the safe zone.
a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
b.prime = (critter.speed^2-a.prime^2)^.5

# This corrects the movement to ensure that the critters are flying at each other rather than away from each other.
movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))
between = function(xy1,xy2,point) (point>xy1&pointxy2&pointmovement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

# Set the new xypos
xypos1 = xypos+movement

points(xypos1, col="red")

# ------------------------------------------------------
# Let's turn this into an animation.

library(animation)

# loopnum = 100; ncritters=40; inertia = .5; show.grid=T; ani.pause=F; plot.fixed=F; plot.centered=F; brownian = F; arrow = T
flocking <- ani.pause="F," arrow="T)" brownian="F," function="" inertia=".5," loopnum="100," ncritters="40," p="" plot.centered="F," plot.fixed="F," show.grid="T,">
  # Generate xy initial positions.
  # xypos will hold the current critter position while
  # xypos0 will hold the position of the critters the previous time.
  xypos = xypos0 = matrix(runif(ncritters*2),ncol=2)-.5

  movement0 = 0

  # Loop though all of the loops.
  for (i in 1:loopnum) {

  # This specifies the range to be graphed.
  if (plot.fixed) rangex=rangey = -.5:.5
  if (!plot.fixed) {
    rangex = c(min(xypos[,1]), max(xypos[,1]))
    rangey = c(min(xypos[,2]), max(xypos[,2]))
  }

  #  This handles the grid size when
  if (plot.centered&!plot.fixed) {
    rangex=c(-max(abs(xypos[,1])), max(abs(xypos[,1])))
    rangey=c(-max(abs(xypos[,2])), max(abs(xypos[,2])))
  }

  # This centers the plot at the middle (0,0) if the plot width is also set to be fixed.
  if (plot.centered&plot.fixed) rangex=rangex-mean(rangex)
  if (plot.centered&plot.fixed) rangey=rangey-mean(rangey)

  # Draw critters
  plot(xypos, main="Swarming Animation", xlab="X", ylab="Y", axes=F, ylim=rangey, xlim=rangex, type="p")
  # Draw arrows.  The start of the arrows is the previous periods location.
  if (arrow&i>1) arrows(x0=xypos0[,1],y0=xypos0[,2],x1=xypos[,1],y1=xypos[,2], length = .1)

  # Show the grid in the background.
  if (show.grid) {
    abline(v=seq(-10,10,.1))
    abline(h=seq(-10,10,.1))
  }

  # Show the origin
  text(0,0, "(0,0)")

  # Calculate each critters distance from each other
  distances = NULL
  for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

  # Drop those within the safe zone.
  distances[abs(distances)
  # This selects the critter closest to the selected critter.
  closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]

#   distances[as.apply(!apply(distances, 1, is.na),1,sum)==0,]=0

  # As done above
  ab = xypos-xypos[closest,]

  a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
  b.prime = (critter.speed^2-a.prime^2)^.5

  movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))

  between = function(xy1,xy2,point) (point>xy1&pointxy2&point
  movement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

  movement[is.na(movement)]=0

  movement0 = movement0*inertia + movement

  # This fancy dodad allows half of the change in movement to be due to random variation.
  if (brownian) movement0=movement0+matrix(rnorm(ncritters*2),ncol=2)*critter.speed/2

  # Set the previous round's xy position to be equal to the current round's.
  xypos0 = xypos

  # Update the current round's.
  xypos = xypos+movement0

  # This is only used in the event that the animate package is in use.
  if (ani.pause) ani.pause()
  }

}

# This generates a GIF animation demonstrating smoothly how these GIFs can be incorper
  ani.options(ani.width=400, ani.height=400, interval=.1)

# You must have imagemagick installed for this to work.
  saveGIF(flocking(300,100,.999, ani.pause=T), movie.name = "Swarming.gif", replace=T)

# Here are two different graphs generated by the previous command(though the one on the bottom uses 200 frames while the one on the top uses 300)

# Let's see how this works.
flocking()
flocking(400,100,.99)
flocking(400,100,.99, plot.fixed=T)

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.