API tools faq
paste
Advertisement
michaelbrules

Untitled

michaelbrules
Sep 16th, 2019
126
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.12 KB | None | 0 0
raw download clone embed print report
  1. from scipy import *
  2. from pylab import *
  3.  
  4. #class containing vector operations. Inherited by ODEsolver
  5. class vecMethods:
  6. def __init__(self):
  7. pass
  8.  
  9. def vecAddition(self, a, b, *args):
  10. if len(a) == len(b):
  11. c=[i+j for i,j in zip(a,b)]
  12. return c
  13. elif len(args) > 0:
  14. partialVec=a[args[1]]
  15. c=[i+j for i,j in zip(partialVec, b)]
  16. return c
  17. else:
  18. raise Exception("Vector length does not match!")
  19.  
  20. def vecUpdate(self, a, b, *args):
  21. if len(args) > 0:
  22. a[args[1]]=self.vecAddition(a[args[1]], b)
  23. return a
  24. else:
  25. return self.vecAddition(a, b)
  26.  
  27. def scalarMulti(self, s, v):
  28. return [s*i for i in v]
  29.  
  30. def absoluteVal(self, a, b):
  31. if len(a) == len(b):
  32. sumR=0
  33. for i in range(len(a)):
  34. sumR+=(a[i]-b[i])**2
  35. return sqrt(sumR)
  36.  
  37. else:
  38. raise Exception("Vector length does not match!")
  39.  
  40.  
  41. #class for solving ODE:s. Inherits vecMethods.
  42. class ODEsolver(vecMethods):
  43. def __init__(self, tMax, y0, numberTimeSteps, a, b, m): #y0 is a vector (list) containing the initial values
  44. self.y0=y0
  45. self.numTimeSteps=numberTimeSteps
  46. self.h=2*tMax/self.numTimeSteps
  47. self.a=a
  48. self.b=b
  49. self.m=m
  50. self.t=[i*self.h for i in range(self.numTimeSteps)] #discretize time
  51.  
  52. def solve(self, fcn, eulerOrRK): #eulerOrRK is a string determining whether to use Euler's or the Runge-Kutta method.
  53. y=[self.y0] #this list will contain the solution to the ODE in the form of y_n:s.
  54. if eulerOrRK == 'e':
  55. for i in range(self.numTimeSteps): #generate the y_n:s using Euler's method
  56. y.append(self.recursionEuler(self.t[i], y[i], fcn))
  57. return y
  58. else:
  59. for i in range(self.numTimeSteps): #generate the y_n:s using the R-K method
  60. y.append(self.recursionRK(self.t[i], y[i], fcn))
  61. return y
  62.  
  63. def solveCoupled(self, fcn, eulerOrRK, extra): #send [j, n, extra] to fcn via recursion
  64. N=len(self.y0)
  65. y=[self.y0]
  66. if eulerOrRK == 'e':
  67. for i in range(self.numTimeSteps):
  68. tmp=[]
  69. for j in range(N):
  70. tmp.append(self.recursionEuler(self.t[i], y[i], fcn, N, j, extra))
  71. y.append(tmp)
  72. return y
  73. else:
  74. for i in range(self.numTimeSteps):
  75. tmp=[]
  76. for j in range(N):
  77. tmp.append(self.recursionRK(self.t[i], y[i], fcn, N, j, extra))
  78. y.append(tmp)
  79. return y
  80.  
  81. def recursionEuler(self, tN, yN, fcn, *args):
  82. f=fcn(tN, yN, *args)
  83. hf=self.scalarMulti(self.h, f)
  84. return self.vecAddition(yN, hf, *args)
  85.  
  86. def recursionRK(self, tN, yN, fcn, *args):
  87. f1=fcn(tN, yN, *args)
  88. k1=self.scalarMulti(self.h, f1)
  89. k2=self.scalarMulti(self.h, fcn(tN+self.m*self.h, self.vecUpdate(yN, self.scalarMulti(self.m, k1), *args), *args))
  90. return self.vecAddition(yN, self.vecAddition(self.scalarMulti(self.a, k1), self.scalarMulti(self.b, k2), *args), *args)
  91.  
  92.  
  93.  
  94.  
  95. def SIS(t, y): #f(t, x) for the sis model
  96. alpha=0.5
  97. beta=1.2
  98. return [-beta*y[0]*y[1]/y[2]+alpha*y[1], beta*y[0]*y[1]/y[2]-alpha*y[1], 0]
  99.  
  100. def SIR(t, y): #f(t, x) for the sir model
  101. alpha=0.05
  102. beta=1.2
  103. gamma=0.02
  104. return [-beta*y[0]*y[1]/y[3]-gamma*y[0], beta*y[0]*y[1]/y[3]-alpha*y[1], gamma*y[0]+alpha*y[1], 0]
  105.  
  106. def SIRwithPopulations(t, y, *args): #args should be [number of populations, current population, omegaList]
  107. alpha=0.06
  108. beta=0.88
  109. gamma=0
  110. n=args[1]
  111.  
  112. iterationList=[i for i in range(args[0])]
  113. del iterationList[args[1]]
  114.  
  115. travelTermsS=0
  116. travelTermsI=0
  117. travelTermsR=0
  118. for m in iterationList:
  119. travelTermsS+=args[2][m][n]*y[m][0]-args[2][n][m]*y[n][0]
  120. travelTermsI+=args[2][m][n]*y[m][1]-args[2][n][m]*y[n][1]
  121. travelTermsR+=args[2][m][n]*y[m][2]-args[2][n][m]*y[n][2]
  122.  
  123. sN=-beta*y[n][0]*y[n][1]/y[n][3]-gamma*y[n][0]+travelTermsS
  124. iN=beta*y[n][0]*y[n][1]/y[n][3]-alpha*y[n][1]+travelTermsI
  125. rN=gamma*y[n][0]+alpha*y[n][1]+travelTermsR
  126. nN=0
  127. return [sN, iN, rN, nN]
  128.  
  129.  
  130.  
  131.  
  132. p_population=2059484
  133. b_population=2462637
  134. m_population=4963349
  135.  
  136. bp_totalpassengers=976100
  137. pm_totalpassengers=2057800
  138. bm_totalpassengers=355700
  139.  
  140.  
  141. wbtop=(bp_totalpassengers*.5)/b_population
  142. wptob=(bp_totalpassengers*.5)/p_population
  143.  
  144. wmtop=(pm_totalpassengers*.5)/m_population
  145. wptom=(pm_totalpassengers*.5)/p_population
  146.  
  147. wbtom=(bm_totalpassengers*.5)/b_population
  148. wmtob=(bm_totalpassengers*.5)/m_population
  149. vecOperators=vecMethods()
  150. popSizes=[m_population, p_population, b_population] #the population sizes (Melbourne, Perth, Brisbane)
  151.  
  152.  
  153. omega=[[0, wmtop, wmtob], [wptom, 0,wptob],[wbtom, wbtop, 0]]
  154. """
  155. for i in range(len(popSizes)): #calculate the omega_nm
  156. tmp=[]
  157. for j in range(len(popSizes)):
  158. if i == j:
  159. tmp.append(0)
  160. else:
  161. denom=vecOperators.absoluteVal(positionVectors[i], positionVectors[j])
  162. tmp.append(omega0*popSizes[i]**xi*popSizes[j]**(xi-1)/(denom**(2+mu)))
  163. omega.append(tmp)
  164. """
  165.  
  166. inst1=ODEsolver(50, [[4784607, 1, 0, 4784608],[2020138, 0, 0, 2020138], [2379724, 0, 0, 2379724]], 10000, .5, .5, 1) #initialize instance
  167. #Infection starts with 1 person in Melbourne
  168. solution2Pop=inst1.solveCoupled(SIRwithPopulations, 'r', omega)
  169.  
  170. pop1Sus=[]
  171. pop1Inf=[]
  172. pop1Res=[]
  173. pop2Sus=[]
  174. pop2Inf=[]
  175. pop2Res=[]
  176. pop3Sus=[]
  177. pop3Inf=[]
  178. pop3Res=[]
  179.  
  180. for i in range(len(solution2Pop)-1):
  181. pop1Sus.append(solution2Pop[i][0][0])
  182. pop1Inf.append(solution2Pop[i][0][1])
  183. pop1Res.append(solution2Pop[i][0][2])
  184. pop2Sus.append(solution2Pop[i][1][0])
  185. pop2Inf.append(solution2Pop[i][1][1])
  186. pop2Res.append(solution2Pop[i][1][2])
  187. pop3Sus.append(solution2Pop[i][2][0])
  188. pop3Inf.append(solution2Pop[i][2][1])
  189. pop3Res.append(solution2Pop[i][2][2])
  190.  
  191. plt.scatter(inst1.t, pop1Sus, s=3, c="b", marker="^", label='Susceptible')
  192. plt.scatter(inst1.t, pop1Inf, s=3, c="r", marker="s", label='Infected')
  193. plt.scatter(inst1.t, pop1Res, s=3, c="g", marker="D", label='Resistant')
  194. plt.legend()
  195. plt.xlabel('Time, years')
  196. plt.ylabel('Number of individuals')
  197. plt.title('Melbourne')
  198. plt.show()
  199.  
  200. plt.scatter(inst1.t, pop2Sus, s=3, c="b", marker="^", label='Susceptible')
  201. plt.scatter(inst1.t, pop2Inf, s=3, c="r", marker="s", label='Infected')
  202. plt.scatter(inst1.t, pop2Res, s=3, c="g", marker="D", label='Resistant')
  203. plt.legend()
  204. plt.xlabel('Time, years')
  205. plt.ylabel('Number of individuals')
  206. plt.title('Perth')
  207. plt.show()
  208.  
  209. plt.scatter(inst1.t, pop3Sus, s=3, c="b", marker="^", label='Susceptible')
  210. plt.scatter(inst1.t, pop3Inf, s=3, c="r", marker="s", label='Infected')
  211. plt.scatter(inst1.t, pop3Res, s=3, c="g", marker="D", label='Resistant')
  212. plt.legend()
  213. plt.xlabel('Time, years')
  214. plt.ylabel('Number of individuals')
  215. plt.title('Brisbane')
  216. plt.show()
  217.  
  218. print("Max people infected in Melbourne at one time" , str(max(pop1Inf)) , "at time" , str((pop1Inf.index(max(pop1Inf))/200)) , "years" )
  219. print("Max people infected in Perth at one time" , str(max(pop2Inf)) , "at time" , str((pop2Inf.index(max(pop2Inf))/200)) , "years" )
  220. print("Max people infected in Brisbane at one time" , str(max(pop3Inf)) , "at time" , str((pop3Inf.index(max(pop3Inf))/200)) , "years" )
  221.  
  222. for i in pop2Inf:
  223. if i>1:
  224. Perth1st=((i/200)*365)
  225. break
  226.  
  227. for i in pop3Inf:
  228. if i>1:
  229. Brisbane1st=((i/200)*365)
  230. break
  231.  
  232.  
  233. print("Time of first infection in Perth, ", str(Perth1st), "days")
  234. print("Time of first infection in Brisbane, ", str(Perth1st), "days")
  235.  
  236. Enm=(3600978.8582057036-3602620.7706316267)/3 # errors at time max I
  237. Enb=(1786964.2076815274-1493251.606199126)/3
  238. Enb=(1786148.5211670923-1786964.2076815274)/3
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!