In 1978 the British Medical Journal reported on an outbreak of influenza at a British boarding school. There were `1000` students. The outbreak began with one infected student.

We want to study the spread of the disease through the population of this school. The total population may be divided into three: 
the infected `(i)`, those who have recovered `(r)`, and
those who are still susceptible `(s)` to get the disease.

We will study the disease on a period of `tm` days. One model of propagation uses 3 differential equations:

```
(1) s'(t) = -b * s(t) * i(t)
(2) i'(t) =  b * s(t) * i(t) - a * i(t)
(3) r'(t) =  a * i(t)
```
where `s(t), i(t), r(t)` are the susceptible, infected, recovered at time `t` and 
`s'(t), i'(t), r'(t)` the corresponding derivatives. 
`b` and `a` are constants:
`b` is representing a number of contacts which can spread the disease and
`a` is a fraction of the infected that will recover. 

We can transform equations `(1), (2), (3)` in finite differences
    (https://en.wikipedia.org/wiki/Finite_difference_method#Example:_ordinary_differential_equation)
    (http://www.codewars.com/kata/56347fcfd086de8f11000014)

```
(I)    S[k+1] = S[k] - dt * b * S[k] * I[k]
(II)   I[k+1] = I[k] + dt * (b * S[k] * I[k] - a * I[k])
(III)  R[k+1] = R[k] + dt * I[k] *a
```

The interval `[0, tm]` will be divided in `n` small intervals of length
`dt = tm/n`.
Initial conditions here could be : `S0 = 999, I0 = 1, R0 = 0`
Whatever S0 and I0, R0 (number of recovered at time 0) is always 0.

The function `epidemic` will return the maximum number of infected
as an *integer* (truncate to integer the result of max(I)).

# Example:

```
tm = 14 ;n = 336 ;s0 = 996 ;i0 = 2 ;b = 0.00206 ;a = 0.41
epidemic(tm, n, s0, i0, b, a) --> 483
```
    
# Notes: 

- You will pass the tests if
`abs(actual - expected) <= 1`

- Keeping track of the values of susceptible, infected and recovered you can plot the solutions of the 3 differential equations. See an example below on the plot.

![alternative text](http://i.imgur.com/xB6VSqzm.png)

Source: codewars | Difficulty: EASY
Tags: Fundamentals

Starter code: def epidemic(tm, n, s0, i0, b, a):
	

Solutions:
def epidemic(tm, n, s, i, b, a):

	def f(s, i, r):
		dt = tm / n
		for t in range(n):
			(s, i, r) = (s - dt * b * s * i, i + dt * (b * s * i - a * i), r + dt * i * a)
			yield i
	return int(max(f(s, i, 0)))


def epidemic(tm, n, s0, i0, b, a):
	timing = 0
	dt = tm / n
	r0 = 0
	i_max = i0
	while timing < tm:
		timing += dt
		s1 = s0 - dt * b * s0 * i0
		i0 = i0 + dt * (b * s0 * i0 - a * i0)
		r0 = r0 + dt * r0 * a
		s0 = s1
		i_max = i0 if i0 > i_max else i_max
	return i_max


def epidemic(tm, n, s0, i0, b, a):

	def evolve():
		(s, i, r, dt) = (s0, i0, 0, tm / n)
		while 1:
			(s, i, r) = (s - dt * b * s * i, i + dt * (b * s * i - a * i), r + dt * i * a)
			yield i
	for i1 in evolve():
		if i1 < i0:
			return i0
		i0 = i1


def epidemic(tm, n, s0, i0, b, a):
	dt = tm / float(n)
	t = [0] * (n + 1)
	s = [0] * (n + 1)
	i = [0] * (n + 1)
	r = [0] * (n + 1)
	s[0] = s0
	i[0] = i0
	r[0] = 0
	t[0] = 0
	for k in range(n):
		t[k + 1] = (k + 1) * dt
		s[k + 1] = s[k] - dt * b * s[k] * i[k]
		i[k + 1] = i[k] + dt * (b * s[k] * i[k] - a * i[k])
		r[k + 1] = r[k] + dt * i[k] * a
	return int(max(i))


def epidemic(tm, n, s0, i0, b, a):
	S = {0: s0}
	I = {0: i0}
	dt = tm / n
	for k in range(n):
		S[k + 1] = S[k] - tm / n * b * S[k] * I[k]
		I[k + 1] = I[k] + tm / n * (b * S[k] * I[k] - a * I[k])
	return int(I[max(I, key=lambda key: I[key])])


epidemic = lambda t, n, s, i, b, a: ([0 for (d, r, m, S, I, R, z) in [(float(t) / n, 0, i, -1, -1, -1, __import__('itertools'))]], [[0 for (S, I, R, s, i, r, m) in [(lambda t: t + (max(m, t[4]),))((s, i, r, s - d * b * s * i, i + d * (b * s - a) * i, r + d * i * a))]] for _ in z.takewhile(lambda x: I < i, z.repeat(0))], m)[-1]


def epidemic(tm, n, s0, i0, b, a):
	(dt, imax) = (tm / n, 0)
	(s, i, r) = (s0, i0, 0)
	for _ in range(n):
		(s, i, r) = (s - dt * b * s * i, i + dt * (b * s * i - a * i), r + dt * i * a)
		imax = max(int(i), imax)
	return imax


def epidemic(tm, n, s0, i0, b, a):
	from math import floor
	(S, I) = ([s0], [i0])
	dt = tm / n
	for _ in range(n):
		s = S[-1] - dt * b * S[-1] * I[-1]
		i = I[-1] + dt * (b * S[-1] * I[-1] - a * I[-1])
		S.append(s)
		I.append(i)
	return int(max(I))