Coverage for aftercovid/models/covid_sird_cst.py: 100%

40 statements  

« prev     ^ index     » next       coverage.py v7.1.0, created at 2024-03-29 03:09 +0100

1# coding: utf-8 

2""" 

3Implementation of a model for epidemics propagation. 

4""" 

5import numpy.random 

6from sklearn.linear_model import ElasticNet 

7from ._base_sir import BaseSIR 

8from .covid_sird import CovidSIRD 

9 

10 

11class CovidSIRDc(BaseSIR): 

12 """ 

13 Inspiration `Modelling some COVID-19 data 

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

15 This model considers that observed data are not the 

16 true ones. 

17 

18 .. math:: 

19 

20 \\begin{array}{rcl} 

21 S &=& S_{obs} (1 + a)\\\\ 

22 I &=& I_{obs} (1 + b)\\\\ 

23 R &=& R_{obs} (1 + c)\\\\ 

24 D &=& D_{obs} 

25 \\end{array} 

26 

27 Where :math:`S`, :math:`I`, :math:`R` are 

28 :math:`R_{obs}` are observed. 

29 hidden, only :math:`S_{obs}`, :math:`I_{obs}`, 

30 The following equality should be verified: 

31 :math:`S + I + R + D = N = S_{obs} + I_{obs} + R_{obs} + D_{obs}`. 

32 We also get from the previous equations: 

33 

34 .. math:: 

35 

36 \\begin{array}{rcl} 

37 dS &=& dS_{obs} (1 + a) = - \\beta \\frac{IS}{N} = 

38 - \\beta \\frac{I_{obs}S_{obs}}{N}(1+a)(1+b) \\\\ 

39 \\Longrightarrow dS_{obs} &=& - \\beta \\frac{I_{obs}S_{obs}}{N}(1+b) 

40 \\end{array} 

41 

42 And also: 

43 

44 .. math:: 

45 

46 \\begin{array}{rcl} 

47 dD &=& dD_{obs} = \\nu I = \\nu I_{obs} (1+b) 

48 \\end{array} 

49 

50 And as well: 

51 

52 .. math:: 

53 

54 \\begin{array}{rcl} 

55 dR &=& dR_{obs} (1 + c) = \\mu I = \\mu (1 + b) I_{obs} \\\\ 

56 \\Longrightarrow dR_{obs} &=& - \\nu I_{obs} \\frac{1+b}{1+c} 

57 \\end{array} 

58 

59 And finally: 

60 

61 .. math:: 

62 

63 \\begin{array}{rcl} 

64 dI &=& dI_{obs} (1 + b) = -dR - dS - dD = 

65 - \\mu \\frac{1 + b}{1+ c} I_{obs} - \\nu (1+b) I_{obs} - 

66 - \\beta I_{obs}\\frac{S_{obs}}{N} (1 + a)(1 + b) 

67 \\\\ 

68 \\Longrightarrow dI_{obs} &=& - \\nu I_{obs} - \\mu I_{obs} 

69 - \\beta I_{obs}\\frac{S_{obs}}{N} (1 + a) 

70 \\end{array} 

71 

72 This model should still verify: 

73 

74 .. math:: 

75 

76 \\begin{array}{rcl} 

77 S_{obs} + I_{obs} + R_{obs} + D_{obs} &=& N = S + I + R + D \\\\ 

78 &=& S_{obs}(1+a) + I_{obs}(1+b) + R_{obs}(1+c) + D_{obs} 

79 \\end{array} 

80 

81 That gives :math:`a S_{obs} + b I_{obs} + c R_{obs} = 0`. 

82 

83 .. runpython:: 

84 :showcode: 

85 :rst: 

86 

87 from aftercovid.models import CovidSIRDc 

88 

89 model = CovidSIRDc() 

90 print(model.to_rst()) 

91 

92 .. exref:: 

93 :title: SIRDc simulation and plotting 

94 

95 .. plot:: 

96 

97 from pandas import DataFrame 

98 import matplotlib.pyplot as plt 

99 from aftercovid.models import CovidSIRDc 

100 

101 model = CovidSIRDc() 

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

103 df = DataFrame(sims) 

104 print(df.head()) 

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

106 ax.set_xlabel("jours") 

107 ax.set_ylabel("population") 

108 r0 = model.R0() 

109 ax.set_title("Simulation SIRDc\\nR0=%f" % r0) 

110 

111 plt.show() 

112 

113 Visual representation: 

114 

115 .. gdot:: 

116 :script: 

117 

118 from aftercovid.models import CovidSIRDc 

119 model = CovidSIRDc() 

120 print(model.to_dot()) 

121 

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

123 common to SIRx models. This model is not really working 

124 better than :class:`CovidSIRD <aftercovid.covid_sird.CovidSIRD>`. 

125 """ 

126 

127 P0 = [ 

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

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

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

131 ('a', -1.5025135094805093e-08, 

132 'paramètre gérant les informations cachées (S)'), 

133 ('b', 1e-5, 'paramètre gérant les informations cachées (I)'), 

134 ('c', 1e-5, 'paramètre gérant les informations cachées (R)'), 

135 ] 

136 

137 Q0 = [ 

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

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

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

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

142 ] 

143 

144 C0 = [ 

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

146 ] 

147 

148 # c = b (nu + mu) / (mu - b * nu) 

149 # (1 + b) / (1 + c) = 1 - nu * b / mu 

150 eq = { 

151 'S': '-beta * (1 + b) * I * S / N', 

152 'I': ('beta * (1 + a) * I * S / N' 

153 '- nu * I - mu * I'), 

154 'R': 'mu * (1 + b) * I / (1 + c)', 

155 'D': 'nu * (1 + b) * I'} 

156 

157 def __init__(self): 

158 BaseSIR.__init__( 

159 self, 

160 p=CovidSIRDc.P0.copy(), 

161 q=CovidSIRDc.Q0.copy(), 

162 c=CovidSIRDc.C0.copy(), 

163 eq=CovidSIRDc.eq.copy()) 

164 

165 def R0(self, t=0): 

166 '''Returns R0 coefficient. 

167 See :meth:`CovidSIRD.R0 <aftercovid.models.CovidSIRD.R0>`''' 

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

169 

170 def correctness(self, X=None): 

171 ''' 

172 Returns :math:`a S_{obs} + b I_{obs} + c R_{obs} = 0`. 

173 

174 :param X: None to use inner quantities 

175 :return: a number 

176 ''' 

177 if X is None: 

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

179 return (X[:, 0] * self['a'] + X[:, 1] * self['b'] + 

180 X[:, 2] * self['c']) / self['N'] 

181 

182 def update_abc(self, X=None, update=True, alpha=1.0, l1_ratio=0.5): 

183 ''' 

184 Updates coefficients *a*, *b*, *c* so that method 

185 :meth:`correctness <aftercovid.models.CovidSIRDc.correctness>` 

186 returns 0. It uses `ElasticNet 

187 <https://scikit-learn.org/stable/modules/generated/ 

188 sklearn.linear_model.ElasticNet.html>`_. 

189 

190 :param X: None to use inner quantities 

191 :param update: True to update to the coefficients 

192 or False to just return the results 

193 :param alpha: see ElasticNet 

194 :param l1_ratio: see ElasticNet 

195 :return: dictionary 

196 ''' 

197 if X is None: 

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

199 X = X / self['N'] 

200 cst = - self.correctness(X) 

201 model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, fit_intercept=False) 

202 model.fit(X, cst) 

203 coef = model.coef_ 

204 res = dict(a=self['a'] + coef[0], b=self['b'] + coef[1], 

205 c=self['c'] + coef[2]) 

206 if update: 

207 self.update(**res) 

208 return res 

209 

210 def rnd(self): 

211 ''' 

212 Draws random parameters. 

213 Not perfect. 

214 ''' 

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

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

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

218 self['a'] = numpy.random.rand() * 1e-5 

219 self['b'] = numpy.random.rand() * 1e-5 

220 self['c'] = numpy.random.rand() * 1e-5 

221 

222 @staticmethod 

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

224 """ 

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

226 

227 :param epsilon: amplitude 

228 :return: new X 

229 """ 

230 return CovidSIRD.add_noise(X, epsilon=epsilon)