Overview

A simple script that recreates the min/max bifurcation diagrams from Hastings and Powell 1991.

Library Functions

Two useful functions are defined in the module bif.py.

   1 import numpy
   2 
   3 def window(data, size):
   4     """A generator that returns the moving window of length
   5     `size` over the `data`
   6 
   7     """
   8     for start in range(len(data) - (size - 1)):
   9         yield data[start:(start + size)]
  10 
  11 
  12 def min_max(data, tol=1e-14):
  13     """Return a list of the local min/max found
  14     in a `data` series, given the relative tolerance `tol`
  15 
  16     """
  17     maxes = []
  18     mins = []
  19     for first, second, third in window(data, size=3):
  20         if first < second and third < second:
  21             maxes.append(second)
  22         elif first > second and third > second:
  23             mins.append(second)
  24         elif abs(first - second) < tol and abs(second - third) < tol:
  25             # an equilibrium is both the maximum and minimum
  26             maxes.append(second)
  27             mins.append(second)
  28 
  29     return {'max': numpy.asarray(maxes),
  30             'min': numpy.asarray(mins)}

The Model

For speed the model is defined in a fortran file and compiled into a library for use from python. Using this method gives a 100 fold increase in speed. The file uses Fortran 90, which makes using f2py especially easy. The file is named hastings.f90.

module model
    implicit none

    real(8), save :: a1, a2, b1, b2, d1, d2

contains 

    subroutine fweb(y, t, yprime)
        real(8), dimension(3), intent(in) :: y
        real(8), intent(in) :: t
        real(8), dimension(3), intent(out) :: yprime

        yprime(1) = y(1)*(1.0d0 - y(1)) - a1*y(1)*y(2)/(1.0d0 + b1*y(1)) 
        yprime(2) = a1*y(1)*y(2)/(1.0d0 + b1*y(1)) - a2*y(2)*y(3)/(1.0d0 + b2*y(2)) - d1*y(2)
        yprime(3) = a2*y(2)*y(3)/(1.0d0 + b2*y(2)) - d2*y(3)
    end subroutine fweb

end module model

Which is compiled (using the gfortran compiler) with the command: f2py -c -m hastings hastings.f90 --fcompiler=gnu95

The Script

   1 import numpy
   2 from scipy.integrate import odeint
   3 import bif
   4 
   5 import hastings
   6 
   7 # setup the food web parameters
   8 hastings.model.a1 = 5.0
   9 hastings.model.a2 = 0.1
  10 hastings.model.b2 = 2.0
  11 hastings.model.d1 = 0.4
  12 hastings.model.d2 = 0.01
  13 
  14 # setup the ode solver parameters
  15 t = numpy.arange(10000)
  16 y0 = [0.8, 0.2, 10.0]
  17 
  18 def print_max(data, maxfile):
  19     for a_max in data['max']:
  20         print >> maxfile, hastings.model.b1, a_max
  21 
  22 x_maxfile = open('x_maxfile.dat', 'w')
  23 y_maxfile = open('y_maxfile.dat', 'w')
  24 z_maxfile = open('z_maxfile.dat', 'w')
  25 for i, hastings.model.b1 in enumerate(numpy.linspace(2.0, 6.2, 420)):
  26     print i, hastings.model.b1
  27     y = odeint(hastings.model.fweb, y0, t)
  28 
  29     # use the last 'stationary' solution as an intial guess for the
  30     # next run. This both speeds up the computations, as well as helps
  31     # make sure that solver doesn't need to do too much work.
  32     y0 = y[-1, :]
  33 
  34     x_minmax = bif.min_max(y[5000:, 0])
  35     y_minmax = bif.min_max(y[5000:, 1])
  36     z_minmax = bif.min_max(y[5000:, 2])
  37 
  38     print_max(x_minmax, x_maxfile)
  39     print_max(y_minmax, y_maxfile)
  40     print_max(z_minmax, z_maxfile)

Cookbook/Theoretical Ecology/Hastings and Powell (last edited 2011-03-09 06:37:28 by Warren Weckesser)