Population Modeling in Python

By | 2012/10/08

One of the courses I enjoyed most in my Ph.D. program was taught by Prof. Kirk Winemiller on population dynamics. There are various collections of models in various languages out there, and multi-model population dynamic applications. But I still think that there is some utility to rolling my own. Since 2009, I’ve gotten more into Python programming, so I thought that I would take a popular class of population dynamic models and produce a Python module to instantiate them.

A long-time standard method in population modeling is the Leslie matrix. This technique applies when one has data about the age structure of a population and produces estimates going forward by using matrix multiplication to go from the population numbers, fecundity, and survivorship numbers to get the estimate of the population in each age class at the next time step.

A similar method is the Lefkovitch approach. This is still based upon matrix operations, but the underlying data involves stages rather than age structure. This sort of model is often used to capture more complex life histories than are tracked in a Leslie matrix model.

The similarities make it straightforward to incorporate both approaches into one supporting Python class.

The following Python module defines the LMatrix class. The dependencies are the Numpy module and the interval module. I used “pip install interval” to get the interval module on my machine. If you run this module in standalone mode, it runs a test of the LMatrix model with a web-accessible example of a Leslie matrix and of a Lefkovitch matrix.

  1. """
  2. popdyn.py
  3.  
  4. Trying out population dynamics in Python.
  5. Wesley R. Elsberry
  6.  
  7. """
  8.  
  9. class LMatrix:
  10.     """
  11.    LMatrix
  12.  
  13.    A support class for Leslie and Lefkovitch matrix use for
  14.    population dynamics.
  15.  
  16.    This is a generic class that allows for an arbitrary number of age
  17.    classes or stages.
  18.    """
  19.  
  20.     def __init__(self,stAges):
  21.         import numpy as num
  22.         import numpy.matlib as M
  23.         from numpy.matlib import rand,zeros,ones,empty,eye
  24.         import interval
  25.  
  26.         """
  27.        In either Leslie age-structured or Lefkovitch stage-
  28.        structured population modeling, the central feature
  29.        is a special matrix representing both fecundity of
  30.        ages/stages and survivorship in each age/stage.
  31.  
  32.        The Leslie age-structured matrix is slightly simpler,
  33.        since each iteration moves the population forward
  34.        by a time step equal to the difference between the
  35.        age classes.
  36.  
  37.        The Lefkovitch stage-structured matrix,
  38.        on the other hand, may have unequal times spent in
  39.        each stage, and thus other elements of the matrix
  40.        represent the fraction of individuals that continue
  41.        to remain in the stage per time step of the model.
  42.        Those lie on the main diagonal.
  43.  
  44.        The matrix in either case is an N-by-N matrix, where
  45.        N is the number of ages or stages (stAges parameter).
  46.        Because most values in the matrix are zero, we'll
  47.        start with that.
  48.        """
  49.  
  50.         self.stAges = stAges  # Keep track of how many age/stage classes there are
  51.         self.m = zeros((self.stAges,self.stAges))
  52.         self.step = 0  # We are at the beginning
  53.         self.popvec = None
  54.         self.survival = None
  55.         self.recurrence = None
  56.         self.fecundity = None
  57.  
  58.     def LM_AddFecundity(self,fvector):
  59.         """
  60.        Method to set fecundity values for an LMatrix.
  61.  
  62.        This is done by setting the first row of the
  63.        matrix to the values in the vector.
  64.  
  65.        A mismatch between the length of the vector and
  66.        the width of the matrix leaves both unchanged.
  67.        """
  68.         if (fvector.shape[0] == self.stAges):
  69.             # Just replace the row
  70.             self.m[0] = fvector
  71.             # Save it in the object
  72.             self.fecundity = fvector
  73.         else:
  74.             print "Mismatch in size: %s vs. %s" % (self.stAges - 1,fvector.shape[0])
  75.  
  76.     def LM_AddSurvival(self,survival):
  77.         """
  78.        Add the values for survival that shift population members
  79.        from one age/stage to the next.
  80.        The values come in as the "survival" vector, a Numpy array.
  81.        They replace values in the m matrix in the diagonal from
  82.        [1,0] to [N-1,N-2].
  83.        """
  84.         if (survival.shape[0] == (self.stAges - 1)):
  85.             for ii in range(1,self.stAges):
  86.                 self.m[ii,ii-1] = survival[ii-1]
  87.             # Save it in the object
  88.             self.survival = survival
  89.         else:
  90.             print "Mismatch in size: %s vs. %s" % (self.stAges - 1,survival.shape[0])
  91.  
  92.     def LM_AddRecurrence(self,recur):
  93.         """
  94.        Add the values for survival of organisms remaining in the same
  95.        stage. This is for stage-structured population models only.
  96.        The input is as the vector recur, and its values replace those
  97.        in the m matrix along the main diagonal from [1,1] to [N-1,N-1].
  98.        """
  99.         if (recur.shape[0] == (self.stAges - 1)):
  100.             for ii in range(1,self.stAges):
  101.                 self.m[ii,ii] = recur[ii-1]
  102.             # Save it in the object
  103.             self.recurrence = recur
  104.         else:
  105.             print "Mismatch in size: %s vs. %s" % (self.stAges - 1,recur.shape[0])
  106.  
  107.     def LM_SetOneRelation(self,fromState,toState, value):
  108.         """
  109.        Method to set a relation that does not fall on the survival
  110.        diagonal or the recurrence diagonal. This is useful for more
  111.        complex stage-structured population modeling where organisms
  112.        from one stage may graduate to multiple other stages at defined
  113.        rates.
  114.        """
  115.         iv = interval.Interval.between(0,self.stAges-1)
  116.         if ((fromState in iv) and (toState in iv)):
  117.             print self.m
  118.             self.m[toState,fromState] = value
  119.             print self.m
  120.    
  121.     def LM_SetPopulation(self,popvector):
  122.         """
  123.        Another central feature of these models is that the size
  124.        of the population is kept in a 1xN column vector. For the
  125.        implementation here, the actual representation is as a
  126.        Numpy array, which has no column vector as such. This will
  127.        be handled in the actual stepping method.
  128.        """
  129.         if (popvector.shape[0] == (self.stAges)):
  130.             self.popvec = popvector
  131.         else:
  132.             print "Mismatch in size: %s vs. %s" % (self.stAges,popvector.shape[0])
  133.  
  134.     def LM_StepForward(self):
  135.         """
  136.        Do the matrix multiplication to obtain the new population
  137.        vector. Retain the previous population vector.
  138.  
  139.        Handle turning population vector into a column vector for the
  140.        multiplication.
  141.        """
  142.         # Convert the population array to a Numpy matrix and transpose it
  143.         # to get the column vector we need. Multiply the L* matrix by
  144.         # the column vector, resulting in a new column vector with the
  145.         # population at the next step.
  146.         nextpopvec = num.mat(self.m) * num.mat(self.popvec).T
  147.  
  148.         # Save the old population vector
  149.         self.lastpopvec = self.popvec
  150.  
  151.         # Replace the population vector with the new one, which means
  152.         # transposing it and converting to Numpy array type
  153.         self.popvec = num.array(nextpopvec.T)
  154.  
  155.         # Track the number of steps taken
  156.         self.step += 1
  157.  
  158.     def LM_TotalPopulation(self):
  159.         """
  160.        Return the total population size. Sums the "popvec" vector.
  161.        """
  162.         if (None != self.popvec):
  163.  
  164.             # Population vector as array multiplied by column vector of 1s is a sum
  165.             t = num.mat(self.popvec) * ones(self.stAges).T
  166.  
  167.             return t[0,0]
  168.         else:
  169.             return 0.0
  170.  
  171.  
  172. if __name__ == "__main__":
  173.     """
  174.     Generic initialization suggested at
  175.     http://www.scipy.org/NumPy_for_Matlab_Users
  176.    """
  177.     # Make all numpy available via shorter 'num' prefix
  178.     import numpy as num
  179.     # Make all matlib functions accessible at the top level via M.func()
  180.     import numpy.matlib as M
  181.     # Make some matlib functions accessible directly at the top level via, e.g. rand(3,3)
  182.     from numpy.matlib import rand,zeros,ones,empty,eye
  183.     # Define a Hermitian function
  184.     def hermitian(A, **kwargs):
  185.         return num.transpose(A,**kwargs).conj()
  186.  
  187.     # Make some shorcuts for transpose,hermitian:
  188.     #    num.transpose(A) --> T(A)
  189.     #    hermitian(A) --> H(A)
  190.     T = num.transpose
  191.     H = hermitian
  192.  
  193.     import interval
  194.  
  195.        
  196.     # Check it against an existing example data set
  197.     # http://www.cnr.uidaho.edu/wlf448/Leslie1.htm
  198.  
  199.     ex1 = LMatrix(4)
  200.  
  201.     fex1 = num.array([0.5, 2.4, 1.0, 0.0])
  202.     ex1.LM_AddFecundity(fex1)
  203.  
  204.     sex1 = num.array([0.5, 0.8, 0.5])
  205.     ex1.LM_AddSurvival(sex1)
  206.  
  207.     pex1 = num.array([20, 10, 40, 30])
  208.     ex1.LM_SetPopulation(pex1)
  209.  
  210.     print pex1
  211.     print ex1.m
  212.     ex1.LM_StepForward()
  213.     print ex1.popvec
  214.    
  215.     # It checks out!
  216.  
  217.     # Another example, this time of a stage-structured population
  218.     # http://www.afrc.uamont.edu/whited/Population%20projection%20models.pdf
  219.     ex2 = LMatrix(3)
  220.  
  221.     fex2 = num.array([0.0, 52, 279.5])
  222.     ex2.LM_AddFecundity(fex2)
  223.  
  224.     sex2 = num.array([0.024, 0.08])
  225.     ex2.LM_AddSurvival(sex2)
  226.  
  227.     rex2 = num.array([0.25, 0.43])
  228.     ex2.LM_AddRecurrence(rex2)
  229.  
  230.     pex2 = num.array([70.0,20.0,10.0])
  231.     ex2.LM_SetPopulation(pex2)
  232.  
  233.     print pex2
  234.     print ex2.m
  235.     ex2.LM_StepForward()
  236.     print ex2.popvec
  237.     print ex2.LM_TotalPopulation()
  238.  
  239.     ex2.LM_StepForward()
  240.     print ex2.LM_TotalPopulation()
  241.  
  242.     ex2.LM_StepForward()
  243.     print ex2.LM_TotalPopulation()
  244.  
  245.     for ii in range(22):
  246.         ex2.LM_StepForward()
  247.     print ex2.popvec
  248.    
  249.     # Tests OK!

Output from the standalone run:

  1. [20 10 40 30]
  2. [[ 0.5  2.4  1.   0. ]
  3.  [ 0.5  0.   0.   0. ]
  4.  [ 0.   0.8  0.   0. ]
  5.  [ 0.   0.   0.5  0. ]]
  6. [[ 74.  10.   8.  20.]]
  7. [ 70.  20.  10.]
  8. [[  0.00000000e+00   5.20000000e+01   2.79500000e+02]
  9.  [  2.40000000e-02   2.50000000e-01   0.00000000e+00]
  10.  [  0.00000000e+00   8.00000000e-02   4.30000000e-01]]
  11. [[ 3835.       6.68     5.9 ]]
  12. 3847.58
  13. 2093.1914
  14. 5811.535142
  15. [[ 19837904.89838918    393232.36554185     30519.85368983]]
Be Sociable, Share!
<> 260965 23618 >

6 thoughts on “Population Modeling in Python

  1. Karen S.

    For anyone who is interested, coursera.org is offering a free online introductory course in Python programming.

  2. JoeBuddha

    Never been a Python fan; never could see the point. More of a C# guy, myself. Maybe I’ll do a translation.

  3. Wesley R. Elsberry Post author

    Language wars! We know from the Invariance Theorem that differences in implementing two programs in different languages with the same output amounts to a fixed constant depending on the two languages chosen, so overall if a language gets you to being able to implement any computable function, it comes down to other factors as to which you choose.

    When I want to shave off the most microseconds I can, I’ll dust off the C++ compiler. For most embedded microcontrollers, I’m more comfortable in C. I still use Perl for the odd script on FreeBSD or Linux, though I’m turning to Python more often now even for those. For LAMP web programming, I use PHP like many others. C# is a natural choice for Microsoft shop environments for desktop, web, and cloud projects. For putting together desktop user interfaces, I used to use Delphi mostly, and now have shifted pretty much to C#. I *like* C#, really, but I’d be lying to say that I’d want to use it for everything I do programming for.

    I have been thinking of a C# version of LMatrix, too. I’ve gotten the Aforge.NET library, but haven’t gotten into it enough yet to do the job. The little bit of matrix multiplication needed could be done in roll-your-own fashion easily enough.

    Python’s strength lies in the toolkits provided as modules. With Numpy, Scipy, and Matplotlib, I’ve got most of what I’d ever want from MATLAB, and in a free, open source tool that anybody else can use, too. There’s lots more available. C# does not have that level of support for scientific programming yet, and add-ons tend far more often to be commercial. Python integrates nicely with all sorts of other tools, including GNU R. If you know COM and R, you can interface C#, too, but my experience is that there is a steeper learning curve more commonly in trying to get things done in C#. I think some of the above considerations made Python the default choice of language for the Raspberry Pi project’s educational focus.

  4. JoeBuddha

    When all you have is a hammer… ;)

    I do all of my work in C# at the moment, although I’ve been known to write the occasional C++ or Java ap. I’ve just never been interested in Python-style languages, and I’ve always developed my own tools if I don’t find them out there. I’m usually programming for fun when this happens, and rolling my own is part of it.

  5. Gerdien

    What would be the easiest program for individual-based programming on a laptop?

Comments are closed.