Population Modeling in Python

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.

Update: The code is available from its own Github repository. That has been tested for compatibility in both Python 2.7 (which is itself past end of life now) and Python 3.7. It should be fine in more recent Python 3 versions, too, but add an issue on the repo if you find a problem. The code here has character transliteration issues, so please snag the tested code from Github instead.

<br />
“””<br />
popdyn.py</p>
<p>Trying out population dynamics in Python.<br />
Wesley R. Elsberry</p>
<p>“””</p>
<p>class LMatrix:<br />
    “””<br />
    LMatrix</p>
<p>    A support class for Leslie and Lefkovitch matrix use for<br />
    population dynamics.</p>
<p>    This is a generic class that allows for an arbitrary number of age<br />
    classes or stages.<br />
    “””</p>
<p>    def __init__(self,stAges):<br />
        import numpy as num<br />
        import numpy.matlib as M<br />
        from numpy.matlib import rand,zeros,ones,empty,eye<br />
        import interval</p>
<p>        “””<br />
        In either Leslie age-structured or Lefkovitch stage-<br />
        structured population modeling, the central feature<br />
        is a special matrix representing both fecundity of<br />
        ages/stages and survivorship in each age/stage.</p>
<p>        The Leslie age-structured matrix is slightly simpler,<br />
        since each iteration moves the population forward<br />
        by a time step equal to the difference between the<br />
        age classes.</p>
<p>        The Lefkovitch stage-structured matrix,<br />
        on the other hand, may have unequal times spent in<br />
        each stage, and thus other elements of the matrix<br />
        represent the fraction of individuals that continue<br />
        to remain in the stage per time step of the model.<br />
        Those lie on the main diagonal.</p>
<p>        The matrix in either case is an N-by-N matrix, where<br />
        N is the number of ages or stages (stAges parameter).<br />
        Because most values in the matrix are zero, we’ll<br />
        start with that.<br />
        “””</p>
<p>        self.stAges = stAges  # Keep track of how many age/stage classes there are<br />
        self.m = zeros((self.stAges,self.stAges))<br />
        self.step = 0  # We are at the beginning<br />
        self.popvec = None<br />
        self.survival = None<br />
        self.recurrence = None<br />
        self.fecundity = None</p>
<p>    def LM_AddFecundity(self,fvector):<br />
        “””<br />
        Method to set fecundity values for an LMatrix.</p>
<p>        This is done by setting the first row of the<br />
        matrix to the values in the vector.</p>
<p>        A mismatch between the length of the vector and<br />
        the width of the matrix leaves both unchanged.<br />
        “””<br />
        if (fvector.shape[0] == self.stAges):<br />
            # Just replace the row<br />
            self.m[0] = fvector<br />
            # Save it in the object<br />
            self.fecundity = fvector<br />
        else:<br />
            print “Mismatch in size: %s vs. %s” % (self.stAges – 1,fvector.shape[0])</p>
<p>    def LM_AddSurvival(self,survival):<br />
        “””<br />
        Add the values for survival that shift population members<br />
        from one age/stage to the next.<br />
        The values come in as the “survival” vector, a Numpy array.<br />
        They replace values in the m matrix in the diagonal from<br />
        [1,0] to [N-1,N-2].<br />
        “””<br />
        if (survival.shape[0] == (self.stAges – 1)):<br />
            for ii in range(1,self.stAges):<br />
                self.m[ii,ii-1] = survival[ii-1]<br />
            # Save it in the object<br />
            self.survival = survival<br />
        else:<br />
            print “Mismatch in size: %s vs. %s” % (self.stAges – 1,survival.shape[0])</p>
<p>    def LM_AddRecurrence(self,recur):<br />
        “””<br />
        Add the values for survival of organisms remaining in the same<br />
        stage. This is for stage-structured population models only.<br />
        The input is as the vector recur, and its values replace those<br />
        in the m matrix along the main diagonal from [1,1] to [N-1,N-1].<br />
        “””<br />
        if (recur.shape[0] == (self.stAges – 1)):<br />
            for ii in range(1,self.stAges):<br />
                self.m[ii,ii] = recur[ii-1]<br />
            # Save it in the object<br />
            self.recurrence = recur<br />
        else:<br />
            print “Mismatch in size: %s vs. %s” % (self.stAges – 1,recur.shape[0])</p>
<p>    def LM_SetOneRelation(self,fromState,toState, value):<br />
        “””<br />
        Method to set a relation that does not fall on the survival<br />
        diagonal or the recurrence diagonal. This is useful for more<br />
        complex stage-structured population modeling where organisms<br />
        from one stage may graduate to multiple other stages at defined<br />
        rates.<br />
        “””<br />
        iv = interval.Interval.between(0,self.stAges-1)<br />
        if ((fromState in iv) and (toState in iv)):<br />
            print self.m<br />
            self.m[toState,fromState] = value<br />
            print self.m</p>
<p>    def LM_SetPopulation(self,popvector):<br />
        “””<br />
        Another central feature of these models is that the size<br />
        of the population is kept in a 1xN column vector. For the<br />
        implementation here, the actual representation is as a<br />
        Numpy array, which has no column vector as such. This will<br />
        be handled in the actual stepping method.<br />
        “””<br />
        if (popvector.shape[0] == (self.stAges)):<br />
            self.popvec = popvector<br />
        else:<br />
            print “Mismatch in size: %s vs. %s” % (self.stAges,popvector.shape[0])</p>
<p>    def LM_StepForward(self):<br />
        “””<br />
        Do the matrix multiplication to obtain the new population<br />
        vector. Retain the previous population vector.</p>
<p>        Handle turning population vector into a column vector for the<br />
        multiplication.<br />
        “””<br />
        # Convert the population array to a Numpy matrix and transpose it<br />
        # to get the column vector we need. Multiply the L* matrix by<br />
        # the column vector, resulting in a new column vector with the<br />
        # population at the next step.<br />
        nextpopvec = num.mat(self.m) * num.mat(self.popvec).T </p>
<p>        # Save the old population vector<br />
        self.lastpopvec = self.popvec</p>
<p>        # Replace the population vector with the new one, which means<br />
        # transposing it and converting to Numpy array type<br />
        self.popvec = num.array(nextpopvec.T)</p>
<p>        # Track the number of steps taken<br />
        self.step += 1</p>
<p>    def LM_TotalPopulation(self):<br />
        “””<br />
        Return the total population size. Sums the “popvec” vector.<br />
        “””<br />
        if (None != self.popvec):</p>
<p>            # Population vector as array multiplied by column vector of 1s is a sum<br />
            t = num.mat(self.popvec) * ones(self.stAges).T</p>
<p>            return t[0,0]<br />
        else:<br />
            return 0.0</p>
<p>if __name__ == “__main__”:<br />
    “””<br />
     Generic initialization suggested at<br />
     http://www.scipy.org/NumPy_for_Matlab_Users<br />
    “””<br />
    # Make all numpy available via shorter ‘num’ prefix<br />
    import numpy as num<br />
    # Make all matlib functions accessible at the top level via M.func()<br />
    import numpy.matlib as M<br />
    # Make some matlib functions accessible directly at the top level via, e.g. rand(3,3)<br />
    from numpy.matlib import rand,zeros,ones,empty,eye<br />
    # Define a Hermitian function<br />
    def hermitian(A, **kwargs):<br />
        return num.transpose(A,**kwargs).conj()</p>
<p>    # Make some shorcuts for transpose,hermitian:<br />
    #    num.transpose(A) –> T(A)<br />
    #    hermitian(A) –> H(A)<br />
    T = num.transpose<br />
    H = hermitian</p>
<p>    import interval</p>
<p>    # Check it against an existing example data set<br />
    # http://www.cnr.uidaho.edu/wlf448/Leslie1.htm</p>
<p>    ex1 = LMatrix(4)</p>
<p>    fex1 = num.array([0.5, 2.4, 1.0, 0.0])<br />
    ex1.LM_AddFecundity(fex1)</p>
<p>    sex1 = num.array([0.5, 0.8, 0.5])<br />
    ex1.LM_AddSurvival(sex1)</p>
<p>    pex1 = num.array([20, 10, 40, 30])<br />
    ex1.LM_SetPopulation(pex1)</p>
<p>    print pex1<br />
    print ex1.m<br />
    ex1.LM_StepForward()<br />
    print ex1.popvec</p>
<p>    # It checks out!</p>
<p>    # Another example, this time of a stage-structured population<br />
    # http://www.afrc.uamont.edu/whited/Population%20projection%20models.pdf<br />
    ex2 = LMatrix(3)</p>
<p>    fex2 = num.array([0.0, 52, 279.5])<br />
    ex2.LM_AddFecundity(fex2)</p>
<p>    sex2 = num.array([0.024, 0.08])<br />
    ex2.LM_AddSurvival(sex2)</p>
<p>    rex2 = num.array([0.25, 0.43])<br />
    ex2.LM_AddRecurrence(rex2)</p>
<p>    pex2 = num.array([70.0,20.0,10.0])<br />
    ex2.LM_SetPopulation(pex2)</p>
<p>    print pex2<br />
    print ex2.m<br />
    ex2.LM_StepForward()<br />
    print ex2.popvec<br />
    print ex2.LM_TotalPopulation()</p>
<p>    ex2.LM_StepForward()<br />
    print ex2.LM_TotalPopulation()</p>
<p>    ex2.LM_StepForward()<br />
    print ex2.LM_TotalPopulation()</p>
<p>    for ii in range(22):<br />
        ex2.LM_StepForward()<br />
    print ex2.popvec</p>
<p>    # Tests OK!</p>
<p>

Output from the standalone run:

<br />
[20 10 40 30]<br />
[[ 0.5  2.4  1.   0. ]<br />
 [ 0.5  0.   0.   0. ]<br />
 [ 0.   0.8  0.   0. ]<br />
 [ 0.   0.   0.5  0. ]]<br />
[[ 74.  10.   8.  20.]]<br />
[ 70.  20.  10.]<br />
[[  0.00000000e+00   5.20000000e+01   2.79500000e+02]<br />
 [  2.40000000e-02   2.50000000e-01   0.00000000e+00]<br />
 [  0.00000000e+00   8.00000000e-02   4.30000000e-01]]<br />
[[ 3835.       6.68     5.9 ]]<br />
3847.58<br />
2093.1914<br />
5811.535142<br />
[[ 19837904.89838918    393232.36554185     30519.85368983]]</p>
<p>