# Rate of growth of aphid population

May I know how to speed up this simulation process? It is simulating a sequence with if condition at each step.

The theta here contains two elements lambda and mu, which are both scalars determining the rate of growth of aphid population, where $$\N[t]\$$ is the current population and $$\C[t]\$$ the accumulative population (including those died). At each time step, the population increase/decrease with probability $$\theta[0] * N[i] * dt\$$ and $$\theta[1] * N[i] * C[i] * dt\$$, otherwise remain the same.

import numpy as np

def aphid(theta):
dt = 0.0001
d  = np.array([10, 20])
N  = [28]
C  = [28]
i  = 0
for t in np.arange(0+dt, np.max(d)+dt, dt):
U = np.random.uniform(0, 1, 1)
if U < theta[0] * N[i] * dt:
N.append(N[i] + 1)
C.append(C[i] + 1)
elif U < theta[0] * N[i] * dt + theta[1] * N[i] * C[i] * dt:
N.append(N[i] - 1)
C.append(C[i])
else:
N.append(N[i])
C.append(C[i])
i=i+1

integers = [int(i) for i in np.floor(d/dt)]
out      = np.array(N)[integers]

return out

• Thanks for the comments and I have made the editions. Basically I am generating a sequence with added elements at each time step a random variable Nov 18 at 20:50
• What are typical values of theta? Nov 20 at 13:17

Welcome to Code Review.

I'll come to speed in a moment, but first I wanted to pick out a few other aspects of this code.

First, I'd suggest you'd benefit from being a little bit more verbose. If you're following a standard formula then there could be value to matching variable names to those used in the formula. If that's the case, it's worth having a comment referencing and explaining that formula. In more general code, single letter variables make things hard to read. That's even more true when one variable (theta) holds two meaningfully different values. It's much cleaner to just have separate variables for lambda and mu, or "breeding_probability" and "dying_probability" or whatever they're supposed to be. Although "integers" is not a single letter name, it is also not hugely descriptive of why you want to use it. Meanwhile some of the magic values that you're using could probably do with being pulled out into named parameters. For example instead of a random [28] ininitialising N, it could be [inital_population].

Second, looking at the simulation itself, there are a couple of things which I'm nervous about, and one that I wanted to focus on. To be clear, I am not an entomologist and I recognise that aphid reproduction is weird. I certainly don't know much about what impact individual specimens being pregnant with a cloned granddaughter would have on population growth! However, from a purely mathematical perspective, I'm a bit nervous about the way you're using your random input. It looks to me like you're just adding up probabilities, which is a bit suspect unless they're guaranteed to be very small. In doing that, you're introducing a dependence in the death check on the birth check. Again, for sufficiently small probabilities this may not matter. Consider, however, if the population was so large that theta[0] * N[i] * dt exceeded 1. That would mean that exactly 1 aphid gets born each step and that the aphids become immortal! Clearly in practice that would be an unreasonable number of aphids, but it should still be a prompt that something isn't quite right. It's worth either doing the maths to work out what the right probability distributions actually are, or quantifying the assumptions made by the model and explicitly checking that any inputs don't break those assumptions.

Actually, I'm curious now. While looking at the actual simulation, do you know the reason that the death rate depends on the number of aphids which have ever existed? I can understand it depending on N, but depending on N*C surprises me.

The final comment I should make is about seeding. If you're doing any sort of random simulation for scientific purposes, it is worth using some sort of random seed at the start. The goal is to ensure that you can re-run the exact same simulation if needed, whether that's for a paper or for making it easier to debug. If you want to see the effect of different simulations with the same parameters (e.g. mu and lambda) you'd just use different seed values.

Anyway, all that out of the way I should talk a bit about code speed.

As always, the first rule of optimising code is to measure it. One of my first thoughts was to replace the integers and out lines with out = np.array([N[int(i)] for i in np.floor(d/dt)]). That is, in your current penultimate line you're constructing a numpy array of two hundred thousand elements only to throw away almost all of them. With my suggestion, all that work gets skipped, and skipping work is the most elegant way to get done faster! However, if I'd measured your code (which I haven't, but you should) I'd find that that final array construction takes negligible time compared with the actual simulation. It's better to get a 5% acceleration on something which takes 99% of the time than to make something which takes only 0.1% of the time unimaginably faster. So once you run something like cProfile or (which I'd actually suggest in this case) line_profiler, you'll be better placed to focus your efforts.

Once you've profiled the code and know where the time is being spent, here are a few tips you may want to use to sort them out:

• If a significant amount of time is being spent redoing the same calculations, such as theta[0] * dt, consider pulling that into a variable you can calculate once.
• If a significant amount of time is spent in the various calls to append consider pre-allocating the arrays to their full length. Every time you call append, there is a risk that Python goes "I don't have enough space. I need to copy this entire array to a bigger bit of free memory." By pre-allocating it, you avoid all those copies. Alternatively, consider whether you actually need to store all that data in the first place. (Again, the fastest work is the work you don't do)
• If a significant amount of time is being spent in the np.uniform call, consider whether you can exploit numpy's optimisations for doing a large batch of work in one go and pre-generate all your random numbers rather than getting one at a time.
• If none of these are enough, you'll need to do a bit more maths to avoid a chunk of the work. For example, I would expect (at least if you're in the small probabilities regime I mentioned earlier) that a majority of your time steps there will be no change. As such, you could instead use a geometric distribution to find out how many no-change steps you'll have before the next thing that actually changes. You'd then need a different equation to find out whether that change is a death or a birth, before repeating the time to the next event. To be clear, this would be a significant change to the nature of your simulation even if the results are the same.
• In all of this, keep re-running your profiler as you make the changes to make sure they do actually help.