Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# coding: utf-8 

2""" 

3Implementation of a model for epidemics propagation. 

4""" 

5import numpy.random 

6from ._base_sir import BaseSIR 

7 

8 

9class CovidSIRD(BaseSIR): 

10 """ 

11 Inspiration `Modelling some COVID-19 data 

12 <http://webpopix.org/covidix19.html>`_. 

13 

14 .. runpython:: 

15 :showcode: 

16 :rst: 

17 

18 from aftercovid.models import CovidSIRD 

19 

20 model = CovidSIRD() 

21 print(model.to_rst()) 

22 

23 .. exref:: 

24 :title: SIRD simulation and plotting 

25 

26 .. plot:: 

27 

28 from pandas import DataFrame 

29 import matplotlib.pyplot as plt 

30 from aftercovid.models import CovidSIRD 

31 

32 model = CovidSIRD() 

33 sims = list(model.iterate(60)) 

34 df = DataFrame(sims) 

35 print(df.head()) 

36 ax = df.plot(y=['S', 'I', 'R', 'D'], kind='line') 

37 ax.set_xlabel("jours") 

38 ax.set_ylabel("population") 

39 r0 = model.R0() 

40 ax.set_title("Simulation SIRD\\nR0=%f" % r0) 

41 

42 plt.show() 

43 

44 Visual representation: 

45 

46 .. gdot:: 

47 :script: 

48 

49 from aftercovid.models import CovidSIRD 

50 model = CovidSIRD() 

51 print(model.to_dot()) 

52 

53 See :ref:`l-base-model-sir` to get the methods 

54 common to SIRx models. 

55 """ 

56 

57 P0 = [ 

58 ('beta', 0.5, 'taux de transmission dans la population'), 

59 ('mu', 1 / 14., '1/. : durée moyenne jusque la guérison'), 

60 ('nu', 1 / 21., '1/. : durée moyenne jusqu\'au décès'), 

61 ] 

62 

63 Q0 = [ 

64 ('S', 9990., 'personnes non contaminées'), 

65 ('I', 10., 'nombre de personnes malades ou contaminantes'), 

66 ('R', 0., 'personnes guéries (recovered)'), 

67 ('D', 0., 'personnes décédées'), 

68 ] 

69 

70 C0 = [ 

71 ('N', 10000., 'population'), 

72 ] 

73 

74 eq = { 

75 'S': '- beta * S / N * I', 

76 'I': 'beta * S / N * I - mu * I - nu * I', 

77 'R': 'mu * I', 

78 'D': 'nu * I' 

79 } 

80 

81 def __init__(self): 

82 BaseSIR.__init__( 

83 self, 

84 p=CovidSIRD.P0.copy(), 

85 q=CovidSIRD.Q0.copy(), 

86 c=CovidSIRD.C0.copy(), 

87 eq=CovidSIRD.eq.copy()) 

88 

89 def R0(self, t=0): 

90 ''' 

91 Returns R0 coefficient. 

92 

93 :param t: unused 

94 

95 This coefficients tells how many people one 

96 contagious infects in average. 

97 Let's consider the contagious population :math:`I_t`. 

98 A time *t+1*, there are :math:`\\beta I_t` new cases 

99 and :math:`(\\nu + \\mu)I_t` disappear. So at time 

100 *t+2*, there :math:`\\beta I_t(1 - \\nu - \\mu)` 

101 new cases, at time *t+d+1*, it is 

102 :math:`\\beta I_t(1 - \\nu - \\mu)^d`. 

103 We finally find that: 

104 :math:`R_0=\\beta \\sum_d (1 - \\nu - \\mu)^d` and 

105 :math:`R_0=\\frac{\\beta}{\\nu + \\mu}`. 

106 ''' 

107 return self['beta'] / (self['nu'] + self['mu']) 

108 

109 def correctness(self, X=None): 

110 """ 

111 Unused. 

112 """ 

113 if X is None: 

114 X = self.vect().reshape((1, -1)) 

115 return numpy.zeros(X.shape) 

116 

117 def rnd(self): 

118 ''' 

119 Draws random parameters. 

120 Not perfect. 

121 ''' 

122 self['beta'] = numpy.random.randn(1) * 0.1 + 0.5 

123 self['mu'] = numpy.random.randn(1) * 0.1 + 1. / 14 

124 self['nu'] = numpy.random.randn(1) * 0.1 + 1. / 21 

125 

126 @staticmethod 

127 def add_noise(X, epsilon=1.): 

128 """ 

129 Tries to add reasonable noise to the quantities stored in *X*. 

130 

131 :param epsilon: amplitude 

132 :return: new X 

133 """ 

134 rnd = numpy.random.randn(*X.shape) * epsilon + 1. 

135 rnd = numpy.maximum(rnd, 0) 

136 

137 X2 = X.copy().astype(numpy.float64) 

138 grad = (X2[1:] - X2[:-1]) 

139 grad[:, 3] *= rnd[1:, 3] 

140 grad[:, 2] *= (rnd[1:, 2] - 1) / 5 + 1. 

141 

142 X2[1:] = X2[0, :] + numpy.cumsum(grad, axis=0) 

143 

144 fact = numpy.multiply(numpy.sum(X2, axis=1), 1. / numpy.sum(X, axis=1)) 

145 X2 = numpy.multiply(X2, fact.reshape(X.shape[0], 1)) 

146 delta = numpy.sum(X, axis=1) - numpy.sum(X2, axis=1) 

147 delta /= X.shape[1] 

148 X2 = numpy.add(X2, delta.reshape(-1, 1)) 

149 return X2