disease model (fwd)

Damien Joly wrote: 

> I'm a population ecologist interested in wildlife
> disease, and have used modeling extensively in my
> research.

I actually played with some similar problems for a project in my Numeric 
Analasys II class.  It was just comparing different methods for controlling 
a population of Impalas (the animals, not the cars).
It wasn't based too closely to any real situation,
but it was kinda fun.
I was actually pretty suprised to find a copy lying around,
i thought i had deleted it.
Oh, i don't really have any place to put it up on the web,
but i guess it may be useful for it to exist at least once,
so i'll just include it with my message (sorry).
it's three files total. 

oh, and the last little part that tries to come up with a stable killing 
rate doesn't work.
i think it may be oscilatory around the equilibrium point or something like 
that,
cause when i do the calculations by hand (matrix theory, Jacobian analysis, 
etc.),
i come up with rates around 15%,
unlike the 80% i've gotten with the program.
oh well, i still got an A ^_^ 

 ------------------------------>8---------------------------
#!/usr/local/bin/ruby -w 

#    Copyright (c) 2002, Martin Chase
#
#    This software is really boring, as it's my homework.
#    Feel free to use this however, subject to the
#    artistic license, as per perl's at:
#          http://www.perl.com/lang/misc/Artistic.html
# 

#   Impala modeling assignment
#   Assumptions:  females only live for 11 years, males for 10
#                 the initial population is evenly distributed 

require "Rules.rb" 

module Impala 

### Module constants
$MALE = 1
$FEMALE = 2
$NEITHER = 0
$JUVENILE_MALES = (0..0)
$YOUNG_ADULT_MALES = (1..5)
$ADULT_MALES = (6..10)
$JUVENILE_FEMALES = (0..0)
$YOUNG_ADULT_FEMALES = (1..5)
$ADULT_FEMALES = (6..11) 

class ImpalaPop
 attr_accessor :males, :females, :rules 

 protected
 ################################
 ###  initialize
 ###  Arguments: hash of named parameters:
 ###                 'initial' = initial population, to be spread evenly 
throughout
 ###                 'rules' = the rules set to use (Rules)
 ###                 'males' = male population spread (Array)
 ###                 'females' = female population spread (Array)
 def initialize(*args)
   args = args[0]
   if (args.has_key?('initial'))
	initial = args['initial']
	remainder = initial % 21
	@males = Array.new
	@females = Array.new
	0.upto(9) do |age|
	  @males << (initial / 21)
	  if ((remainder -= 1) >= 0);
	    @males[-1] += 1
	  end
	end
	0.upto(10) do |age|
	  @females << (initial / 21)
	  if ((remainder -= 1) > 0)
	  @females[-1] += 1
	  end
	end
   end 

   if (args.has_key?('rules'))
	@rules = args['rules']
   end 

   if (args.has_key?('males'))
	@males = args['males']
   end 

   if (args.has_key?('females'))
	@females = args['females']
   end
 end 

 public
 ################################
 ###  nextYear
 ###  returns a ImpalaPop object representing the population of next years 
impalas,
 ###    based on the rules provided, or defaults.
 def nextYear (rules = nil)
   next_males = Array.new(1)   #first entry undefined
   next_females = Array.new(1) #first entry undefined
   (0..10).each {|i|
	next_males[i+1] = @males[i] * @rules.male_survival[i] if @males[i]
	next_females[i+1] = @females[i] * @rules.female_survival[i]
   }
   births = expected_births
   if (births * 2 > 3 * adult_males)
	births = 3 / 2 * adult_males
   end
### this algorithm matches pam's, but is wrong.  a male is only
### needed to fertilize a female.
#        fertile_females = young_adult_females + adult_females
#        if (fertile_females > adult_males * 3)
#  	births = adult_males * 3
#        else
#  	births = expected_births
#        end 

   next_males[0] = next_females[0] = births  #gender proportional births 

   return ImpalaPop.new('males' => next_males, 'females' => next_females,
			   'rules' => @rules)
 end 

 ################################
 ###  expected_births
 ###  ignoring males, how many female babies should be born
 def expected_births
   births = 0
   @females.each_index {|i|
	births += @females[i] * @rules.birthrate[i]
   }
   return births
 end 

 ################################
 ###  population
 def population (gender = $NEITHER, range = nil)
   total_pop = 0
   @females.each {|f| total_pop += f}
   @males.each {|m| total_pop += m}
   return total_pop
 end 

 ################################
 ###  gives a string of the population in each section of each gender
 def to_s
   r =("Impala population: %d [Juvenile Males: %d, Young Males: %d, Adult 
Males:" +
	  "%d, Juvenile Females: %d, Young Females: %d, " +
	  "Adult Females: %d]") % [population,
	                           juvenile_males,
	                           young_adult_males,
	                           adult_males,
	                           juvenile_females,
	                           young_adult_females,
	                           adult_females]
   return r
 end 

 ################################
 ###  adult_males
 def adult_males
   dracula = 0
   (6..10).each do |i|
	dracula += @males[i] if @males[i]
   end
   return dracula
 end 

 ################################
 ###  young_adult_males
 def young_adult_males
   dracula = 0
   (1..5).each do |i|
	dracula += @males[i]
   end
   return dracula
 end 

 ################################
 ###  juvenile_males
 def juvenile_males
   return @males[0]
 end 

 ################################
 ###  adult_females
 def adult_females
   dracula = 0
   (6..11).each do |i|
	dracula += @females[i]
   end
   return dracula
 end 

 ################################
 ###  young_adult_females
 def young_adult_females
   dracula = 0
   (1..5).each do |i|
	dracula += @females[i]
   end
   return dracula
 end 

 ################################
 ###  juvenile_females
 def juvenile_females
   return @females[0]
 end 

end 

end 

 ----------------------------->8--------------------------- 

#!/usr/local/bin/ruby -w 

#    Copyright (c) 2002, Martin Chase
#
#    This software is really boring, as it's my homework.
#    Feel free to use this however, subject to the
#    artistic license, as per perl's at:
#          http://www.perl.com/lang/misc/Artistic.html
# 

#   Impala modeling assignment
#   Assumptions:  females only live for 11 years, males for 10
#                 the initial population is evenly distributed
#                 one adult male is needed to fertilized up to 3 females 

class NilClass
def coerce (*args)
 return [0,0] #this is to smooth over math by coercing 'nil' into zero
end
def + (arg)
 return arg
end
end 

module Impala 

class Rules
 attr_accessor :male_survival, :female_survival, :birthrate 

 private
 ##################################
 ###  initialize
 ###  arguments: ([male death rates], [female death rates], [birth rates])
 def initialize (males, females, births)
   @male_survival = males
   @female_survival = females
   @birthrate = births
 end 

 public
 ##################################
 ###  alter_survival(!)?
 ###  arguments:
 ###      gender - one of the module constants $MALE or $FEMALE
 ###      age - the age group to be affected (can be a range)
 ###      degree - the (in|de)crease in survival, as expressed by a decimal
 def alter_survival!(gender, age, degree)
   age = (age.is_a?(Integer)) ? [age] : age.to_a
   gender = case gender
		 when $MALE
		   @male_survival
		 when $FEMALE
		   @female_survival
	       end
   degree = degree / age.length #apply equally to ages
   age.each {|a|
	gender[a] += degree
   }
 end 

 def alter_survival(gender, age, degree)
   (self.dup).alter_survival!(gender, age, degree)
 end 

 def clone
   return Rules.new(@male_survival.dup, @female_survival.dup,
		       @birthrate.dup)
 end
end 

end 

 ------------------------>8---------------------------- 

#!/usr/local/bin/ruby -w 

#    Copyright (c) 2002, Martin Chase
#
#    This software is really boring, as it's my homework.
#    Feel free to use this however, subject to the
#    artistic license, as per perl's at:
#          http://www.perl.com/lang/misc/Artistic.html
# 

#   Impala modeling assignment
#   Assumptions:  females only live for 11 years, males for 10
#                 the initial population is evenly distributed 

require "ImpalaPop"
require "Rules" 

class Array 

################################
###  sum - add the values of the array together (using '+')
def sum
 value = self[0]
 self.each_index { |v| value += self[v] unless (v == 0)}
 return value
end 

################################
###  average - add the values and divide by the length
def average
 self.sum / self.length
end 

################################
###  applyAndSum - apply a method to each element, combining the results
###  argument:
###      symbol - the symbol for the method to apply
###  this is intended for use with things like '+', '*', &c., but
###  may have other uses.
def applyAndSum (symbol)
 value = self[0]
 self.each_index { |v|
   value = value.send(symbol,self[v]) unless (v == 0)}
 return value
end
end 

module Impala 

$rules = Rules.new([.6,.8,.95,1,1,1,1,.75,.34,0,0],
		     [.6,.9,.95,.97,.97,.95,.95,.95,.8,.7,0],
		     [0,.35,.45,.45,.45,.45,.45,.45,.45,.45,.45]) 

###################################
###  calculate_finals - get numbers (final_pop, growth_rates,
###                     avg_growth) for a given rule set and
###                     a given amount of time.
###  returns: [ImpalaPop, growth_rates, avg_growth_rate]
###  arguments:
###    rules - the rule set to use
###    iterations - the # of time frames to iterate through
###    initial - the initial population
def calculate_finals(rules, iterations, initial = 220)
 pop = ImpalaPop.new('rules' => rules,
			'initial' => initial)
 growth_rates = Array.new 

 i = 0
 i.upto(iterations) do
   last = pop.population
   pop = pop.nextYear
   growth_rates << (pop.population / last)
 end 

 return [pop, growth_rates, growth_rates.average]
end 

###################################
###  print_results - in the style of the day
###  args:
###    title - name it, baby!
###    rate - how fast can you burn
###    male, female, juvenile, y_adult, adult - rates per group
def print_results(title, rate, male, female, juvenile,
		    y_adult, adult)
 puts "\n#{title}\n\n"
 puts "Same initial population.\n"
 rules = $rules.clone
 rules.alter_survival!($MALE, $JUVENILE_MALES,
			  rate*male*juvenile)
 rules.alter_survival!($MALE, $YOUNG_ADULT_MALES,
			  rate*male*y_adult)
 rules.alter_survival!($MALE, $ADULT_MALES,
			  rate*male*adult)
 rules.alter_survival!($FEMALE, $JUVENILE_FEMALES,
			  rate*female*juvenile)
 rules.alter_survival!($FEMALE, $YOUNG_ADULT_FEMALES,
			  rate*female*y_adult)
 rules.alter_survival!($FEMALE, $ADULT_FEMALES,
			  rate*female*adult)
 iterations = 15
 (pop, growth, avg) = calculate_finals(rules, iterations)
 puts "\nAfter #{iterations} years:\n" + pop.to_s 

 puts "\nGrowth rates:\n" + growth.inspect
 puts "Average: " + avg.to_s
 puts "\n\n" 

end 

###################################
###  stabalize - find the population management rate that
###              stabalizes the population growth.  prints.
###  arguments:
###    title - label me.
###    iterations - the number of years to test at
###    tolerance - the allowable error in the growth rate (a
###                decimal < .01, for any useful info).
###    step_size - the starting step size to alter the rate.
###    starting_rate - the initial rate to start testing at.
###    male, ..., adult - the specific rates to apply to the
###                       population.
def stabalize (title, iterations, tolerance, step_size, starting_rate,
				 male, female, juvenile, y_adult, adult)
 puts "\nEquilibrium rate: #{title}\n"
 puts "\nWithin #{tolerance} error, the equilibrium solution is:\n" 

 rate = starting_rate - step_size
 sign_change = false
 (a, growths, c) = calculate_finals($rules, 15)
 last = growths[-1]
 growth = nil 

 until ((tolerance > (1 - last).abs) || (last == growth)) do
   break if rate > 0  #prevent badness 

   step_size /= 2 if sign_change
   rate = rate - ((1-last>0) ? -step_size : step_size )
   rules = $rules.clone
   rules.alter_survival!($MALE, $JUVENILE_MALES,
			    rate*male*juvenile)
   rules.alter_survival!($MALE, $YOUNG_ADULT_MALES,
			    rate*male*y_adult)
   rules.alter_survival!($MALE, $ADULT_MALES,
			    rate*male*adult)
   rules.alter_survival!($FEMALE, $JUVENILE_FEMALES,
			    rate*female*juvenile)
   rules.alter_survival!($FEMALE, $YOUNG_ADULT_FEMALES,
			    rate*female*y_adult)
   rules.alter_survival!($FEMALE, $ADULT_FEMALES,
			    rate*female*adult)
   (pop, growth, avg) = calculate_finals(rules, iterations)
   growth = growth[-1]
   sign_change = ( ((last > 1) && (growth < 1)) ||
                   ((last < 1) && (growth > 1)) ) ? true : false
#      puts rate.to_s + " -- " + growth.to_s + ":" + last.to_s +
#           " -- " + step_size.to_s + "." + sign_change.to_s
   last, growth = growth, last
 end 

 puts rate.to_s + "\n\n"
end 

end 

if __FILE__ == $0 

	include Impala
	
	puts "IMPALA " * 4 + "PROJECT!!!!\n\n"
	
	
	### UNREGULATED GROWTH
	puts "\nUnregulated growths\n"
	pop = ImpalaPop.new('rules' => $rules,
						'initial' => 220)
	puts "\nInitial population:\n" + pop.to_s
	
	iterations = 15
	(pop, growth_rates, avg) = calculate_finals($rules, iterations)
	
	puts "\nAfter #{iterations} years:\n" + pop.to_s
	
	puts "\nGrowth rates:\n" + growth_rates.inspect
	puts "\nAverage: " + avg.to_s
	puts "\n\n"
	
	### PREDATION - LOW
	title = "Predation - Low"
	rate     = -.06 # six percent predation rate
	male     = .5
	female   = .5
	juvenile = .45
	y_adult  = .2
	adult    = .35
	print_results(title, rate, male, female, juvenile, y_adult, adult)
	
	### PREDATION - HIGH
	title = "Predation - High"
	rate     = -.16 # sixteen percent predation rate
	male     = .5
	female   = .5
	juvenile = .45
	y_adult  = .2
	adult    = .35
	print_results(title, rate, male, female, juvenile, y_adult, adult)
	
	### TROPHY HUNTING - LOW
	title = "Trophy Hunting - Low"
	rate     = -.06 # six percent predation rate
	male     = .7
	female   = .3
	juvenile = .025
	y_adult  = .025
	adult    = .95
	print_results(title, rate, male, female, juvenile, y_adult, adult)
	
	### TROPHY HUNTING - HIGH
	title = "Trophy Hunting - High"
	rate     = -.16 # sixteen percent predation rate
	male     = .7
	female   = .3
	juvenile = .025
	y_adult  = .025
	adult    = .95
	print_results(title, rate, male, female, juvenile, y_adult, adult)
	
	### GAME RANCHING - LOW
	title = "Game Ranching - Low"
	rate     = -.06 # six percent predation rate
	male     = .7
	female   = .3
	juvenile = .05
	y_adult  = .75
	adult    = .2
	print_results(title, rate, male, female, juvenile, y_adult, adult)
	
	### GAME RANCHING - HIGH
	title = "Game Ranching - High"
	rate     = -.16 # six percent predation rate
	male     = .7
	female   = .3
	juvenile = .05
	y_adult  = .75
	adult    = .2
	print_results(title, rate, male, female, juvenile, y_adult, adult)
	
	
	###############################
	### E Q U I L I B R I U M S ###
	###############################
	
	iterations = 15
	tolerance = .001
	step_size = .01
	starting_rate = -.01
	
	###  PREDATION
	title    = "Predation"
	male     = .5
	female   = .5
	juvenile = .45
	y_adult  = .2
	adult    = .35
	stabalize(title, iterations,  tolerance, step_size, starting_rate,
			  male, female, juvenile, y_adult, adult)
	
	###  TROPHY HUNTING
	title    = "Trophy Hunting"
	male     = .7
	female   = .3
	juvenile = .025
	y_adult  = .025
	adult    = .95
	stabalize(title, 50,  tolerance, step_size, starting_rate,
			  male, female, juvenile, y_adult, adult)
	#:!: use 50 above because 15 wasn't producing feasable results
	
	###  GAME RANCHING
	title    = "Game Ranching"
	male     = .7
	female   = .3
	juvenile = .05
	y_adult  = .75
	adult    = .2
	stabalize(title, iterations,  tolerance, step_size, starting_rate,
			  male, female, juvenile, y_adult, adult)
	
end