FAQ and troubleshooting

General questions

How does matpopmod compare to other libraries?

There are several R packages to study matrix population models, among which popbio, popdemo, Rage, and Rcompadre. Matpopmod implements most of the functionalities provided by those packages, as well as a few additional ones.

The main difference between matpopmod and the packages listed above is that matpopmod is meant to be used with Python rather than with R. However, it is also possible to use matpopmod with R, as explained below.

Another difference is that matpopmod was designed from the start to work with large datasets, where it is not possible to inspect each model manually. As a result, it is “stricter” regarding the enforcement of mathematical assumptions: for instance, if a formula is valid for primitive projection matrices only, then matpopmod will refuse to apply it to an imprimitive matrix and, when asked to do so, will issue a message detailing the problem. This is meant to help you catch potential errors in your code or in your data.

Matpopmod also has a well-documented API with a modular structure that make it easy to add functionalities. So if you find that something is missing, please visit our GitLab Repository to see how it can can be added to the next release.

Finally, another software with a scope similar to that of matpopmod is ULM. Matpopmod owes a lot to ULM and to its author Stéphane Legendre. The main specificity of ULM is that it is aimed at users with little to no experience with programming languages, and comes with a fully fledged graphical user interface. However, this comes at the price of less flexibility in the manipulation of models. Moreover, being written in a declining programming language, the future of ULM is uncertain.

Can I use matpopmod with R?

Yes! After you have installed matpopmod for Python, you can use reticulate to call it from the R console:

> library("reticulate")
> matpopmod <- import("matpopmod")
> m = matpopmod$examples$orcinus_orca
> m$A
       [,1]   [,2]   [,3]   [,4]
[1,] 0.0000 0.0043 0.1132 0.0000
[2,] 0.9775 0.9111 0.0000 0.0000
[3,] 0.0000 0.0736 0.9534 0.0000
[4,] 0.0000 0.0000 0.0452 0.9804
> m$lmbd
[1] 1.025441
> m$elasticities
           [,1]        [,2]       [,3] [,4]
[1,] 0.00000000 0.001513103 0.04069515    0
[2,] 0.04220825 0.336325827 0.00000000    0
[3,] 0.00000000 0.040695151 0.53856251    0
[4,] 0.00000000 0.000000000 0.00000000    0

NumPy array are automatically converted to R vector/matrices, so the interface should feel very familiar to R users. Note that, conversely, Python objects have to be appropriately converted. Thus, when creating a model,

MPM(A = [[0, 2], [1/2, 3/4]],
    metadata = {"ModelName": "Test model"})

will have to be replaced by

MPM(A = list(list(0, 2), list(1/2, 3/4)),
    metadata = list(ModelName = "Test model"))

although in general it will be easier to supply A using an existing matrix or a string:

> my_matrix
     [,1] [,2]
[1,]  0.0 2.00
[2,]  0.5 0.75
> matpopmod$MPM(A = my_matrix)
MPM with projection matrix:
[0.   2.  ]
[0.5  0.75]

Reticulate also makes it possible to use matpopmod with R Markdown and to run Python scripts and get their output in the R console. See the project homepage for more information.

How can I cite matpopmod?

You can cite matpopmod as

Bienvenu, F. and Doulcier, G. (2021). MatPopMod, a Python library for matrix population models. Zenodo. DOI: 10.5281/zenodo.5557426.

Note that this DOI will always point to the latest version of matpopmod. If you want to refer to a specific version of matpopmod, you can do so as follows:

Bienvenu, F. and Doulcier, G. (2021). MatPopMod, a Python library for matrix population models (version 0.1.0). Zenodo. DOI: 10.5281/zenodo.5557427.

The complete list of versions with the corresponding release dates and DOIs can be found on Zenodo.

For LaTeX users, we suggest using the following BibTeX entry

@article{matpopmod,
  title={{M}at{P}op{M}od, a {P}ython library for matrix population models},
  author={Bienvenu, Fran{\c{c}}ois and Doulcier, Guilhem},
  doi={10.5281/zenodo.5557426},
  journal={Zenodo},
  year={2021},
}

In addition to this, please consider citing the scientific publication for which matpopmod was developed (in prep).

I found a bug, what can I do?

Please get in touch with us!

The most convenient for us would be if you could open an issue on our issue tracker and describe the problem there in as much detail as possible. If you can, try to include the following information:

  • What operating system are you using?

  • What version of Python and NumPy are you using?

  • What are the minimal steps to reproduce the problem?

  • What output were you expecting and what did you get?

If you do not have a GitLab account and do not wish to create one, you can also just send us an email.

Can I contribute to matpopmod? Request features?

Sure! Although of course we cannot guarantee that feature requests will get implemented. See the previous section for how to contact us.

Troubleshooting

I get a different result using another library

This section is about when the numerical values returned by matpopmod differ from those returned by other libraries. Situations where matpopmod and other libraries disagree about whether a descriptor is defined are discussed in the next section.

Make sure that the discrepancy is not simply due to a difference of convention. In particular, note that the vector or reproductive values v returned by matpopmod is scaled so that \(\mathbf{vw} = \sum_i v_i w_i = 1\), where w is the stable stage distribution, whereas popbio scales it so that \(v_1 = 1\) and ULM so that \(\sum_i v_i = 1\).

That said, the numerical values returned by matpopmod can sometimes differ significantly from those returned by other libraries / software. In every such case that we could identify, the correct values were those returned by matpopmod.

Fortunately, the situations in which such discrepancies arise are rare and almost always correspond to matrices that are either non-conventional (e.g, reducible, periodic, with an unusual ordering of the classes…) or exceptional in some sense (such as \(\lambda = 1\)). For instance, consider the following Leslie model:

\[\begin{split}\begin{pmatrix} 0 & 1 & 1 \\ \frac{1}{3} & 0 & 0 \\ 0 & \frac{2}{3} & \frac{2}{3} \\ \end{pmatrix}\end{split}\]

Here, the asymptotic growth rate is \(\lambda = 1\), so one has to be careful that the \(R_0\) generation time cannot be computed with the formula \(T_{R_0} = \log R_0 / \log \lambda\). Matpopmod says that \(T_{R_0} = 4\), which is the correct value [Elln18].

>>> m = MPM(S = [[0, 0, 0], [1/3, 0, 0], [0, 2/3, 2/3]],
...         F = [[0, 1, 1], [0, 0, 0], [0, 0, 0]])
>>> m.T_R0
3.9999999999999996

However, popbio gives a different numerical value:

> A = rbind(c(0, 1, 1), c(1/3, 0, 0), c(0, 2/3, 2/3))
> generation.time(A)
[1] 0.6666667

To see why this value is absurd and why \(T_{R_0} = 4\), consider what happens when we modify the entries of A infinitesimally:

> A_bis = rbind(c(0, 1, 1 + 1e-12), c(1/3, 0, 0), c(0, 2/3, 2/3))
> generation.time(A_bis)
[1] 4.002667
> A_ter = rbind(c(0, 1, 1 - 1e-12), c(1/3, 0, 0), c(0, 2/3, 2/3))
> generation.time(A_bis)
[1] 4.002667

Note that here the value of \(T_{R_0}\) returned by popbio was “wrong enough” to clearly indicate that there was a problem. However, this is not always going to be the case:

> A_quat = rbind(c(0, 1, 1 + 1e-15), c(1/3, 0, 0), c(0, 2/3, 2/3))
> generation.time(A_quat)
[1] 3

(note that here matpopmod still returns the correct value \(T_{R_0} = 4\)).

One situation in which the numerical value computed by matpopmod will often differ from that computed by another library is the calculation of the period of second-order oscillations by ULM. Here discrepancies are due to the fact that ULM uses an incorrect mathematical definition (see second_order_period for matpopmod’s definition). For instance, for the model dipsacus_sylvestris ULM gives a period of 6.276, whereas the period computed by matpopmod is 2.936:

>>> dipsacus = matpopmod.examples.dipsacus_sylvestris
>>> dipsacus.second_order_period
2.935524088519757

Simulating the trajectory confirms that the period of the second-order oscillations of largest amplitude is indeed 2.926. The following code will display the figure below.

dipsacus = mpm.examples.dipsacus_sylvestris
n0 = [0, 0, 0, 0, 0, 1] # initial population vector
traj = dipsacus.trajectory(n0, t_max=25)
mpm.plot.trajectory(
  traj,
  second_order = True,
  rescale = True,
  show_period = True
)
_images/faq-1.png

Why does matpopmod return nan here?

Sometimes, matpopmod will say that some descriptors are not well-defined and return nan, even though other libraries might give a numerical value for those descriptors. This is mainly going to occur with models that are not primitive, and is going to be accompanied by a warning the first time that this happens with a given model. For instance,

>>> from matpopmod.examples import bernardelli_beetle
>>> bernardelli_beetle.A
array([[0.        , 0.        , 6.        ],
       [0.5       , 0.        , 0.        ],
       [0.        , 0.33333333, 0.        ]])
>>> bernardelli_beetle.w # stable stage distribution
UserWarning: A is not quasi-primitive. Most descriptors are
ill-defined. They will be set to NaN.
array([nan, nan, nan])

Yet libraries such as popbio will return a numerical value:

> stable.stage(A)
[1] 0.6 0.3 0.1

This difference is a design choice. Here, mathematically there is no stable stage distribution because the population is going to oscillate indefinitely. What popbio returns is the Perron vector, i.e. the unique positive eigenvector associated to the dominant eigenvalue of an irreducible matrix (see [Meye00], Section 8.3). We could have chosen to return the Perron vector, but our documentation says that w is the stable stage distribution and some of our subsequent calculations rely on this interpretation – so we prefer to return nan. If the Perron vector is needed, it can be obtained with right_eigenvectors[0]:

>>> bernardelli_beetle.right_eigenvectors[0]
array([0.6, 0.3, 0.1])

This design choice is not mere mathematical nitpicking: it is meant to make help the user catch potential errors in their code. This is especially useful when working with large amounts of matrix population models that cannot be inspected individually: in that situation, asking for the stable distribution of a periodic model is likely to mean that we have forgotten that this model could be periodic, or maybe even that there is an error in the model.

Finally, note that there are also situations in which matpopmod will return a numerical value whereas other libraries will return nan or NA. In that case, this means that we are using more general formulas, often based on recent theoretical results. For instance, matpopmod is able to compute the \(R_0\) generation time of the following matrix

>>> m = MPM(S = "0 0; 0.4 0.6", F = "0 1; 0, 0")
>>> m.T_R0
3.4999999999999996

whereas popbio is not:

> generation.time(A)
[1] NaN

The reason why matpopmod is able to compute \(T_{R_0}\) here even though from the formula \(T_{R_0} = \log R_0 / \log \lambda\) it may look like it is not defined when \(\lambda = 1\), is that, as shown by [Elln18], as \(\lambda \to 1\) we have \(T_{R_0} \to T_a\), the mean age of mothers,

How can I suppress those warnings?

As discussed in the previous section, matpopmod will sometimes issue warnings when it thinks something should be brought to the attention of the user. In general, we recommend not turning off warnings. However, if you have a specific reason to do so then you can use:

import warnings

with warnings.catch_warnings():
    warnings.simplefilter("ignore", category = UserWarning)
    # your code here

Another possibility is to turn those warnings into exceptions and use them to control the flow of your program. See the documentation of Python’s warnings module for details.

Why do I get ValueError: assignment destination is read-only?

The following code will, perhaps puzzlingly, throw an error:

>>> from matpopmod.examples import orcinus_orca
>>> orcinus_orca.v
array([ 1.14163164,  1.19762278,  1.79386902, -0.        ])
>>> rescaled_v = orcinus_orca.v
>>> for i in range(len(rescaled_v)):
...     rescaled_v[i] /= orcinus_orca[0]
...
Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
ValueError: assignment destination is read-only

It may seem like our intentions were good: we were merely trying to compute a rescaled vector of reproductive values, with the reproductive value of newborns set to 1. However, this error is perfectly justified. It is in fact a feature that was implemented to limit the risk of errors when using matpopmod.

To understand what is going on, recall that it is not possible to modify the descriptors of a matrix population model in any way:

>>> orcinus_orca.v = numpy.array([0, 0, 0, 0])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: can't set attribute
>>> orcinus_orca.v[0] = 0
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: assignment destination is read-only

This is to ensure that these descriptors are always correct, no matter what programming errors are made when using the library. Now, in Python as in most programming languages, rescaled_v = orcinus_orca.v does not copy orcinus_orca.v, it merely gives it another name. Thus, rescaled_v[i] = x is rigorously equivalent to orcinus_orca.v[i] = x. Had we been allowed to modify rescaled_v, we would also have modified orcinus_orca.v – which may have been a problem if other parts of our code used orcinus_orca.v and we had forgotten about it.

Now, to achieve what we want we can simply use:

rescaled_v = orcinus_orca.v / orcinus_orca.v[0]

which is not only shorter, but also more efficient. The reason why this works is that the operation orcinus_orca.v / x creates a new array, instead of modifying the existing one. If we have to set the entries of rescaled_v one at a time, we can do so as follows:

>>> rescaled_v = orcinus_orca.v.copy()
>>> for i in range(len(rescaled_v)):
...     rescaled_v[i] /= orcinus_orca.v[0]
...
>>> rescaled_v
array([ 1.        ,  1.04904483,  1.57132034, -0.        ])

Why do I get UnexpectedMathError?

If you ever get an UnexpectedMathError, please try to save the model that led to it and contact us: those occur when a mathematical property is violated as a result of roundoff errors in numerical calculations. While a theoretical possibility, we do not expect these errors to occur in practice – and we have never seen one. Thus, either you were incredibly lucky or this could be indicative of a more serious problem.

Why do I get an error with this 1×1 model?

Trying to create a 1×1 model as follows will raise an error:

>>> MPM(A = 1.5)
Traceback (most recent call last):
  ...
matpopmod.utils.IncorrectDims: Matrix A is not 2-dimensional: 0-D.

and so will other commands such as MPM(A = numpy.array(1.5)).

The reason for this is that we need a 2-dimensional object to describe the projection matrix A, such as a list of lists. The following will work:

>>> MPM(A = [[1.5]])
MPM(
  A = [[1.5]],
  metadata = {}
)

And as usual you can also use a string to specify A:

>>> MPM(A = "1.5")
MPM(
  A = [[1.5]],
  metadata = {}
)