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)

Difficulty: introductory

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
        
    # your code


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):
    s = [s0]; i = [i0]; r = [0]
    t=0
    dt = tm/n
    t += dt
    
    for k in range(n):
        s[k] = s[k-1] - dt*b*s[k-1]*i[k-1]
        i[k] = i[k-1] + dt * (b * s[k-1] * i[k-1] - a * i[k-1])
        r[k] = r[k-1] + dt * i[k-1] *a
    
        
        s.append(s[k])
        i.append(i[k])
        r.append(r[k])
    
    
    
    return int(max(i))

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))