1

I am running a research on the impacts of bed nets on mosquito populations and how it affects malaria transmission. EIR tells us the number of infectious bites a person gets in a year. To calculate this, I have to a summation in python which makes use of a recursisve function S(t). The series converges after about 50 iterations but it takes a very long time (about an hour) if I increase the summation limit beyond 20. I might be doing something wrong, I will very much appreciate your help. The code written in python is displayed below.

##IMPORTING PACKAGES
import sympy as sy
import matplotlib.pyplot as plt
import numpy as np
import profile

##Declaring the constants.
f = 1/3       ## frequency of feeding under zero net coverage
tau_1 = 0.69  ## time spent seeking blood under zero net coverage
tau_2 = 2.31  ## time spent resting during feeding cycle
Q_0 = 0.38    ## proportion of blood meals on humans under zero net coverage
mu_M0 = 0.096 ## daily mosquito mortality under zero net coverage
p_1 =0.91     ##probability of a mosquito survivng a blood meal at zero net coverage
p_2 = 0.74    ##probability of surviving resting phase at zero net coverage
phi_LLIN = 0.89   ##proportion of bites taken on humans when in bed
d_LLIN = 0.41     ##probability that a mosquito is killed by LLIN
s_LLIN = 0.03     ##probability that a mosquito feeds successfully with LLIN
r_LLIN = 0.56     ##probability that a mosquito is repelled by LLIN
phi_IRS = 0.97    ##proportion of bites taken on humans while indoors
psi = 0.86

#chi_LLIN = sy.symbols('chi_LLIN')
chi_LLIN=np.linspace(0.0,1.0,101)
w_LLIN =1. - Q_0*chi_LLIN*phi_LLIN*(1.-s_LLIN);w_LLIN

z_LLIN = Q_0 *chi_LLIN*phi_LLIN*r_LLIN; z_LLIN
tau_1chiLLIN = tau_1 /(1.-z_LLIN);tau_1chiLLIN
f_chiLLIN =1./ (tau_1chiLLIN + tau_2);f_chiLLIN
p1_chiLLIN = (p_1*w_LLIN)/(1.-(z_LLIN*p_1));p1_chiLLIN
p_chiLLIN = (p1_chiLLIN*p_2)**f_chiLLIN;p_chiLLIN
mu_MchiLLIN = -np.log(p_chiLLIN);mu_MchiLLIN

##Variable human blood index
Q_LLIN = Q_0 * ((1.+phi_LLIN*chi_LLIN*(s_LLIN-1.))/(1.+Q_0*phi_LLIN*chi_LLIN*(s_LLIN-1.) ))


##DEFINING E
totalx =39 ##limit of the summation
def E1(x):
    E =0.
    for t in xrange(1,x+1):
        E = E+((p_chiLLIN)**t)/f_chiLLIN
    return E


##DEFINING THE SPOROZOITE INFECTION PREVALENCE
n = 11 #the threshold
kappa=0.0297 ##human infectivity to mosquitoes

def S(x):
    if x<=n:    
        return 0
    elif x==n+1:    
        return 0.05
    else:
        return S(x-1) + ((kappa*Q_0*(1-S(x-1)))/f_chiLLIN)

B=0
for x in xrange(1,totalx+1):
    B = B + ((p_chiLLIN)**x)/f_chiLLIN*S(x)


beta_h = (Q_0/f_chiLLIN)*B

EIR = beta_h * E1(totalx)

EIR
user2357112 supports Monica
  • 215,440
  • 22
  • 321
  • 400
Ozymandais
  • 23
  • 5

1 Answers1

4
def S(x):
    if x<=n:    
        return 0
    elif x==n+1:    
        return 0.05
    else:
        return S(x-1) + ((kappa*Q_0*(1-S(x-1)))/f_chiLLIN)
#              ^^^^^^                  ^^^^^^

In S(x), you make 2 recursive calls to compute S(x-1), which make 4 recursive calls to compute S(x-2), which make 8 recursive calls to compute S(x-3), and so on until you hit the base case of the recursion. You don't need to recompute things so much; just compute S(x-1) once and use the value twice:

def S(x):
    if x<=n:    
        return 0
    elif x==n+1:    
        return 0.05
    else:
        s_lower = S(x-1)
        return s_lower + ((kappa*Q_0*(1-s_lower))/f_chiLLIN)

This will still do more recomputation than necessary, though, because every time you call S:

for x in xrange(1,totalx+1):
    B = B + ((p_chiLLIN)**x)/f_chiLLIN*S(x)

it'll have to recompute all values of S down to the base case again. Consider further changes, such as memoizing S.

Community
  • 1
  • 1
user2357112 supports Monica
  • 215,440
  • 22
  • 321
  • 400