% This short project studies a bit about house sparrows in England
% It has been found that house sparrows are declining due to more
% intensive agriculture (D. G. Hole et al, 2002, Widespread local
% house-sparrow extinctions, Nature 418:931).
% These birds follow the matrix studied in class. The paper gives
% rough values of the parameters.

m     = 3.0;  %Roughly the number of females per female
sigma = 0.6;  %Juvenile survivorship to the end of the summer
slow  = 0.35; %Winter survival, low value 
shigh = 0.65; %Winter survival, high value 
% SYNTAX NOTE: The semicolon at the end of the line supresses output

Alow  = [0 m; sigma*slow slow]   %Matrix with low winter survival
Ahigh = [0 m; sigma*shigh shigh] %Matrix with high winter survival
% The juvenile to adult survival multiplies summer and winter survivorships
% SYNTAX NOTE: The square brackets tell matlab you are building a matrix,
% and the semicolon inside the square brackets tells it to start a new row

% Use matlab to find those eigenvalues!
eigval_Alow  = eig(Alow)
eigval_Ahigh = eig(Ahigh)
% The low population is predicted to decline, while the high population
% will grow

% Find the eigenvectors
[eigvec_Alow,Eigval_Alow]   = eig(Alow)    
[eigvec_Ahigh,Eigval_Ahigh] = eig(Ahigh)
% The columns of the matrix evlow are the eigenvectors.  
% Try the command "help eig" to see how to get different outputs from the
% same command

% To find SAD, one associated with the leading eigenvalue
% Notice that the larger eigenvalue appears second in each case 
SAD_Alow  = eigvec_Alow(:,2)/sum(eigvec_Alow(:,2))
SAD_Ahigh = eigvec_Ahigh(:,2)/sum(eigvec_Ahigh(:,2))
% SYNTAX NOTE: eigvec_Alow(:,2) tells matlab to give all rows (the :) in 
% the second column (the 2) of eigvex_Alow. If the leading eigenvalue happened
% to show up first, you'd have to put a 1 there.

% See the population converge to the eigenvector
onead=[0; 1];                                %An initial population with 1 adult
Alow_1gen     = Alow*onead                   %After 1 generation
Alow_1gen_n   = Alow_1gen/sum(Alow_1gen)     %Normalized after 1 generation
Alow_10gen    = Alow^10*onead                %After 10 generations
Alow_10gen_n  = Alow_10gen/sum(Alow_10gen)   %Normalized after 10 generations
Alow_100gen   = Alow^100*onead               %After 100 generations
Alow_100gen_n = Alow_100gen/sum(Alow_100gen) %Normalized after 100 generations

% Compare the value of an adult to a juvenile
onejuv=[1; 0];                               %An initial population with 1 juvenile
adjuv_compare = sum(Alow^100*onead)/sum(Alow^100*onejuv)
% The total population size after 100 generations starting from one adult 
% compared to the total population size after 100 generations starting from 
% one juvenile.   This is the relative reproductive value of the adult.

% Make graphs of the population size over the first 50 generation
allgen_onead=onead;                             %Initialize a matrix to save data
allgen_onejuv=onejuv;                           %Initialize a matrix to save data
for i=1:49,                                     %Start up a loop
  allgen_onead=[allgen_onead Alow^i*onead];     %Save ith value in ith column
  allgen_onejuv=[allgen_onejuv Alow^i*onejuv];  %Save ith value in ith column
end                                             %End of loop
plot(sum(allgen_onead))                         %Plots total population size
% SYNTAX NOTE: sum adds up along columns
plot(1:50,sum(allgen_onead),'b-',1:50,sum(allgen_onejuv),'g--')
% SYNTAX NOTE: The plot command in general likes to read arguments in threes: 
% the first argument (1:50 is shorthand for the numbers from 1 to 50) is the input, 
% the second is the output, and the third is the style (blue line).  
% The next set of three tells it to use the same input, and plot the total
% population size starting from one juvenile using a dashed green line.

% Plot the adult numbers against the juvenile number
plot(allgen_onead(1,:),allgen_onead(2,:),'bo')
% Plots the first row as the horizontal coordinate and the second row as the
% vertical coordinate using blue circles
plot(allgen_onead(1,:),allgen_onead(2,:),'bo',allgen_onejuv(1,:),allgen_onejuv(2,:),'gx')
% Both the adult and juvenile populations converge to the same direction

% The value of adjuv_compare (or the ratio of the heights of the blue and green 
% curves in the graph of total population size) is supposed to be related to 
% the left eigenvector of the matrix, or the eigenvector of the transpose.  
% Let's check
Alow_transpose = Alow';  %The transpose command
[eigvec_Alow_transpose,Eigval_Alow_transpose]   = eig(Alow_transpose)    
% Normalize so the juvenile is worth 1
reproductive_value  = eigvec_Alow_transpose(:,2)/eigvec_Alow_transpose(1,2)
% Adult component of reproductive_value should equal adjuv_compare

% Finally, a bit of sensitivity analysis.  If we could add a wee bit to one
% parameter, which would we choose?
slowup = slow+0.05;               %Increase winter survival
Aup = [0 m; sigma*slowup slowup]; %Matrix with increased, but low, winter survival
slowup_eigval_change = max(eig(Aup))-max(eig(Alow)) %Change in the eigenvalue

sigmaup = sigma+0.05;             %Increase juvenile summer survival
Aup = [0 m; sigmaup*slow slow];   %Matrix with increased juvenile summer survival
sigmaup_eigval_change = max(eig(Aup))-max(eig(Alow))%Change in the eigenvalue

mup = m+0.05;                     %Increase reproduction
Aup = [0 mup; sigma*slow slow];   %Matrix with increased reproduction
mup_eigval_change = max(eig(Aup))-max(eig(Alow))    %Change in the eigenvalue