%% Lab 3 -- Zombie Attack
%% Zombie invasion!
% The plot below is from a simulation of a zombie attack on the MSU campus. The MSU community number around 57,000 total. Each point in the grid represents one individual. Dark red means infected, light red means healthy. As you can see, the MSU community is unlikely to survive past a month under this attack.
%
%
% The following code simulates this zombie attack with the assumptions discussed in the pre-lab.

% Set the total population
total_population = 57600

% Create an array to represent the population
%    The code below creates a matrix that is 1 row by
%    total_population columns, that is filled with zeros
% We will use zero to represent a healthy individual. 
% We change the value to one when the individual is infected. 
infected_status = zeros(1,total_population)

% Select a random person to infect initially
%     randi(x) returns a random integer between 1 and x
% The code below sets the infected status of a randomly chosen
%  person to 1
infected_status(randi(total_population)) = 1

% The total number of infected people can be found by
sum(infected_status)

% To simulate what happens in one day: for each zombie already
%    present, we choose a random person (regardless of whether he 
%    or she is already infected) and change his or her status to 1. 
% For simplicity we allow the zombie to choose to bite himself. 
for i = 1:sum(infected_status)
    infected_status(randi(total_population)) = 1; % Semicolon to suppress output
end

% So in day 2 the number of infected people is
sum(infected_status)

% Now let us run it forward for 10 more days, printing out the number of infected 
% people each day
for j = 1:10   % For loop runs a total of ten times, simulating ten days
    for i = 1:sum(infected_status)  % During each day, each of the zombies gets a bite
        infected_status(randi(total_population)) = 1; 
    end
    sum(infected_status)   % Shows the number of infected people after the zombies feasted
end

% For the first five days, the number of zombies pretty much doubles every day, exhibiting the classical exponential growth type behavior.
% After the 7th or 8th day, however, the population no longer doubles each day, diverging from the exponential growth behavior.
% In fact, if you run the code multiple times, you can get different results from the simulation, but the growth rate always seem to drop (instead of doubling, it increases by something less than 100%). Let's examine this in more detail.
%
% This time, we run the same code, but instead of showing the zombie population at every step, we show the rate of growth compared to the previous day.

total_population = 57600
infected_status = zeros(1,total_population)
infected_status(randi(total_population)) = 1

% Write down on the side how many people are infected on day 1
saved_infected_num = sum(infected_status)   

% Days 2 through 8
for j = 2:8
    for i = 1:sum(infected_status)
        infected_status(randi(total_population)) = 1; 
    end
    % The following line gives the % change from the previous day
    % Remember, rate of change is difference divided by the previous data
    (sum(infected_status)- saved_infected_num) / saved_infected_num * 100
    % Having computed the rate of change, we no longer need the info about
    % the number of people infected at the end of the previous day. But, we
    % need to prepare for tomorrow and write down the number of infected
    % people right now. 
    saved_infected_num = sum(infected_status);
end

% We see that the first 6 or 7 days the rate of growth is 100%, meaning that the population doubles every day.
% Now let's see what happens between days 9 through 14.

% Days 9 through 14
for j = 9:14
    for i = 1:sum(infected_status)
        infected_status(randi(total_population)) = 1; 
    end
    % The following line gives the % change from the previous day
    (sum(infected_status)- saved_infected_num) / saved_infected_num * 100
    saved_infected_num = sum(infected_status);
end

% By the end the rate drops down to around 90%.
% Let's run the code a little bit longer and see what happens in subsequent days

% Days 15 through 22
for j = 15:22
    for i = 1:sum(infected_status)
        infected_status(randi(total_population)) = 1; 
    end
    % The following line gives the % change from the previous day
    (sum(infected_status)- saved_infected_num) / saved_infected_num * 100
    saved_infected_num = sum(infected_status);
end

% We see that the percentage growth rate drops rapidly.
%% Exercise 1
% The following code plots two percentages as a function of the number of days since the "patient zero" appeared. In blue is the percentage grow rate of the zombie community. In red is the percentage of the overall MSU community that is infected.

total_population = 57600
infected_status = zeros(1,total_population)
infected_status(randi(total_population)) = 1

saved_infected_num = sum(infected_status);

% Days 2 through 30
for j = 2:30
    for i = 1:sum(infected_status)
        infected_status(randi(total_population)) = 1; 
    end
    % The following line gives the % change from the previous day
    percent_growth_rate = (sum(infected_status)- saved_infected_num) / saved_infected_num * 100;
    saved_infected_num = sum(infected_status);
    percent_zombified = saved_infected_num / total_population * 100;
    
    % Plot in blue, the point (j,percent_growth_rate), signifying that the
    % zombie population is growing at percent_growth_rate on day #j
    plot(j,percent_growth_rate, 'bo')
    hold on
    % Similarly, in red, plot the percent of zombified individuals. 
    plot(j,percent_zombified, 'ro')
end
hold off

% Discuss with your teammates, referring to
% basic facts we know about zombies
% the graph generated from the code above
% and answer the following two questions.
% Why is the percentage growth rate dropping over time?
% Hypothesize a relation between the "percentage growth rate of the zombie population" with the "percentage of MSU community infected".
%
%% Logistic Equation
% Remember that the differential equation describing exponential growth is
% \frac{dY}{dt} = r \cdot Y
%  where r is the "rate of growth" and is given as a constant. Also recall in the pre-la where we derive the differential equation describing logistic growth as
% \frac{dZ}{dt} = r_0 \cdot \frac{P - Z}{P} \cdot Z
% Where r_0 is the initial rate and P is the number of Spartans overall.
% Comparing this logistic equation with the exponential growth equation, we can imagine the role of the "rate of growth" r in the exponential growth equation to be played by the (variable value) r_0 \cdot \frac{P-Z}{P}. In other words, the "rate of growth" now depends on how many zombies there already are.
% When there is only one individual infected, (P-Z)/P \approx 1, so initially the rate of growth is given by
% r_0 \cdot \frac{P-Z}{P} \approx r_0
% So what is the value of r_0 we should use? The initial growth rate is one that doubles the zombie population every day (see the simulation!)
%
%% Exercise 2
% Find the value of r_0 corresponding to the zombie population doubling every day, using what you know from exponential growth. (Or, more concretely, if Y solves \frac{dY}{dt} = r_0 Y and Y(0) = 1 and Y(1) = 2, find r_0.)
%% Exercise 3
% Using the value of r_0 you found above, and the value of P = 57600 given, solve the differential equation
% \frac{dZ}{dt} = r_0 \cdot \frac{P - Z}{P}\cdot Z
% with initial value such that when t = 1, Z = 1 by following the steps below.
% Step 1
% Remember that r_0 and P are just numbers. They do not depend on t or on Z. So using what you learned from solving differnetial equations in course, you want to move everything depending on Z to one side of the equation, and everything depending on t to the other. Since nothing on the right hand side depends on t, we rewrite as
% \frac{P}{(P-Z)Z}~ dZ = r_0 ~ dt
% Step 2
% The right hand side integrates easily.
% To integrate the left hand side, you may use the identity
% \frac{1}{P-Z} + \frac{1}{Z} = \frac{P}{(P-Z) Z}
% so you end up with
% \int \frac{1}{P-Z} + \frac{1}{Z} ~dZ = \int r_0 ~dt
% Perform this integration and solve for t as a function of Z. Determine the unknown constant of integration using the fact that when t = 1, Z = 1.
% Note
% The algebraic identity that is used in Step 2 is called a "Partial Fraction Decomposition". You will learn later in this course (Calculus II) how to use partial fractions as a general technique for integration.
% Numerical Connection (optional)
% You may want to check your work for exercises 2 and 3 by plotting your solution against the simulation. You can use the following code block for that.)

% Fill in the right hand side of this variable with the value of 
%   r_0 you found
r_0 = 1 % <---- REPLACE BY THE CORRECT VALUE
P = 57600

% Fill in the right hand side of this function definition with
%  the function t(Z) you found in exercise 3.
t = @(Z)  Z  % <----- REPLACE BY THE CORRECT FUNCTION


% You shouldn't need to edit anything below. 

% Plot by simulation
infected_status = zeros(1,P)
infected_status(randi(P)) = 1

plot(1,1,'ro')
hold on
% Days 2 through 30
for j = 2:30
    for i = 1:sum(infected_status)
        infected_status(randi(P)) = 1; 
    end
    plot(j,sum(infected_status), 'ro')
end

% Plot the solution you found
Zs = 1:P
plot(t(Zs),Zs, 'b')
hold off

%
%% Exercise 4
% As the zombie population approaches the total population (or, as Z becomes very close to P), the number of surviving non-zombie individuals (P-Z) decays exponentially. At what rate is this exponential decay? (In other words, (P-Z) can be written as approximately  C e^{-r t}, what is the correct r?) Justify your answer.
%
%% Bonus Material: Spartan Versus Zombies
% Good news, MSU discovered a Zombie cure/vaccine. Once injected, it cures the patient of zombie-ism and gives limited protection to future zombie bites. The effectiveness of the protection decreases by 5% each day, and when the effectiveness drops below 20% of its original value, the person can be infected again. (So the vaccine is effective for just about 1 month.)
% To eradicate the disease, MSU is sending out volunteers to innoculate random members of the community. How many shots must be given per day to wipe out zombies from campus?
% The code below simulates just such a situation. You will find that
% If too few shots are given out per day, this will stem the tide a bit, either by slowing down the rate of infection, or by limiting the maximum number of people infected at any given time. In either of these cases the zombie population will plateau off around a critical number and stay there.
% When the number of zombies can be controlled to be at a small number, we say that the community has gained Herd Immunity. Our zombie vaccine loses effectiveness overtime, so it is more like the flu vaccine; one must constantly get booster shots to control the spread of the disease.
% If enough shots are given out per day, however, the zombie population will rise initially (before the immunization effort has taken full effect), and then slowly drop off. When this happens, we say that we have achieved Eradication through Herd Immunity. This is very difficult to achieve on a world-wide leve. The only documented case of successful eradication of an infectious disease is smallpox, for which no cases have been documented since 1978.

total_population = 57600;
infected_status = zeros(1,total_population);
vaccine_status = zeros(1, total_population);
number_to_jab_per_day = 100;

infected_status(randi(total_population)) = 1;

for j = 2:100
    % Zombies infect people
    for i = 1:sum(infected_status)
        bitten_person = randi(total_population);
        % Check if the person is vaccinated
        if vaccine_status(bitten_person) < 0.2
            infected_status(bitten_person) = 1;
        end 
    end
    % Protection from vaccine age away
    vaccine_status = vaccine_status .* 0.95;
    % Give random people shots
    for i = 1:number_to_jab_per_day
        person_to_jab = randi(total_population);
        infected_status(person_to_jab) = 0;  % Cures the person
        vaccine_status(person_to_jab) = 1;   % Makes the person immune
    end
    % Show the number of zombies
    sum(infected_status)
end

%
%% Appendix
% Just in case you are wondering how the animated gif at the top of the document is created:
% Starting from the simulation code we had before, we just need to add some additional visualization code.

total_population = 57600

infected_status = zeros(1,total_population)

infected_status(randi(total_population)) = 1

% To visualize this on a grid, first convert our array of 
% infected_status into a square matrix using the reshape() 
% command, and we can then visualize it using the spy() command
spy(reshape(infected_status,[240,240]))
% Notice that 57600 = 240^2

 
for i = 1:sum(infected_status)
    infected_status(randi(total_population)) = 1; % Semicolon to suppress output
end

% Visualize Day 2
spy(reshape(infected_status,[240,240]))

% Now let us run it forward for 10 more days, printing out the number of infected 
% people each day
for j = 1:10
    for i = 1:sum(infected_status)
        infected_status(randi(total_population)) = 1; 
    end
end

% Visualize Day 12
spy(reshape(infected_status,[240,240]))