"Hormesis.py -- model facultative delay in reproduction"
pdelay=0.5       # Probability of delaying reproduction under stress
import sys, random, math
random.seed()
print "Hormesis: maternally-inherited delayed reproduction w/ toxin exposure"
# Genotypes
types=["now",'later','maybe']     # This list doesn't grow.
# Initial values.
genotype=types                      # These three lists grow: one element per individual.
age=[0,0,0]                         # Incremented until death, for each individual.
state=["juvenile","juvenile",'juvenile']
births=[0,0,0]                      # Births recorded for each individual.
# Expand lists to 2**n individuals of each genotype to start.  Lists grow with reproduction below.
n=10
for i in range(n):
   genotype=genotype+genotype       # Appending to itself (+ operation on lists) doubles it.
   age=age+age                 
   state=state+state
   births=births+births
live={'now':2**n, 'later':2**n, 'maybe':2**n}     # Initialize population counts, etc. 
dying={'now':0., 'later':0., 'maybe':0.}
fecundity={'now':0., 'later':0., 'maybe':0.}
lifespan={'now':0., 'later':0., 'maybe':0.}
print "year,  stress,  now,  maybe, later, lifen,  lifem,  lifel,  repron,  reprom,  reprol"
# Print header to paste results into spreadsheet.
# Parameters
maturity={'now':2, 'later':3,'maybe':2}
risk={'juvenile':0.3, 'adult':0.2, 'fertile':0.25, 'dead':0}
# Juvenile risk changed below, varying sinusoidally over years.
period=15
cycles=2

year=0
while year < period*cycles :                       # Loop over years.
   year=year+1
   stress=1+math.sin(2*math.pi*(year+8)/period)    # Stress (famine?) varies sinusoidally from 0-2.
   risk['juvenile']= risk['adult']+stress/4        # Stress increases risk to juveniles only.
   for x in types:             # Zero today's stats for each genotype.
       live[x]=0
       dying[x]=0.
       fecundity[x]=0.
       lifespan[x]=0.
 
   for i in range(len(age)):   # Loop over individuals.
      if state[i]<> 'dead' :
          age[i]=age[i]+1                            # Update age unless dead.
          live[genotype[i]]=live[genotype[i]]+1      # Count individuals alive today.
          if age[i]>maturity[genotype[i]]:           # Reproductive maturity?
             state[i]='fertile'
          if genotype[i]=='maybe' and stress>1 and pdelay>random.random() :
              # Genotype "maybe" may delay reproduction when stressed.
              state[i]="adult" 
      if state[i]=="fertile" :                       # Reproducing?
          births[i]=births[i]+1
          genotype.append(genotype[i])
          age.append(0)
          state.append("juvenile")
          births.append(0)
      if risk[state[i]]>random.random() :            # Dying?
          state[i]='dead'
          dying[genotype[i]]=dying[genotype[i]]+1    # Stats for individuals dying this year.
          fecundity[genotype[i]]=fecundity[genotype[i]]+births[i]
          lifespan[genotype[i]]=lifespan[genotype[i]]+age[i]      
 
   # Report results for year
   lifen=lifespan['now']/dying['now']
   lifem=lifespan['maybe']/dying['maybe']
   lifel=lifespan['later']/dying['later']
   repron=fecundity['now']/dying['now']
   reprom=fecundity['maybe']/dying['maybe']
   reprol=fecundity['later']/dying['later']
   print 11*"%8.2f," % (year, risk['juvenile'],live['now'],live['maybe'],live['later'],lifen,lifem,lifel, repron,reprom,reprol)