% mat2valvecs_08  May, 2008
%   calls function file "senstuff.m" 

% reference: Caswell Chapter four.
% Run this program

%Demographic lessons
% eigenvector analysis of an stage-structured matrix
% asymptotic dynamics
% eigenvalues
% eigenvalue spectrum
% real and complex eigenvalues: absolute magnitude
% dominant eigenvalue

%This example is adapted from Caswell's (2001) Chapter 4 (p. 72-79), 
% plus additional investigations of the complex plane 

clear all
close all
%this clears previous workspace

A=[0 0 0 4.665 61.896; 0.675 0.703 0 0 0; 0 0.047 0.657 0 0; 0 0 0.019 0.682 0; 0 0 0 0.061 .8091] %enter a matrix by rows
name=['\it {Caretta caretta}, Crowder et al. 1994']; % creates a string variable

%make a matfile called Caretta that has A in it
save Caretta A

[vec,val] = eig(A)
% finds the eigenvectors and eigenvalues
%  try "help eig" to see more about it

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Practice reading and interpretting the output from eig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a)  There is a list of eigenvalues; there are more eigenvalues than the dominant one!
% Some are real numbers and some are complex numbers.  
% Recall that a complex number is written a+bi, 
% where a is the real part and bi is the imaginary part and i is the square root of -1. 
% Complex numbers occur in conjugate pairs: a+bi and a-bi.  
% Identify the real and the complex eigenvalues.

% b)  There are several eigenvectors.  Each eigenvector is associated with its own eigenvalue.
% Real eigenvalues have real eigenvectors, while complex eigenvalues have complex eigenvectors. 
% Eigenvectors are “directions” and eigenvalues are the “speeds” associated with each direction. 
% Both complex and real numbers affect dynamics. 
% To understand the dynamics over time, one needs to consider how “fast” the system 
%    will move in each direction. 

% c) Write the solution to the projection equation (Caswell 2001,  Equation 4.49, p. 75) for the turtles.  
%  Don't try to find a numerical representation of the c's, just use the symbols for the c's.
%  The rest of the numbers can be found from a consideration of what is the output from the
%     eig function. 
% Give a word-meaning for the projection equation.  
% Save the eigenvalues and eigenvectors of the turtle in a  “*.mat” file, 
%  so you can call them up again if you need them.

save Caretta   %adds the eigenvalues and eigenvectors to the mat file

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Meanwhile, let's focus on the dominant speed and direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Let lambda1 be the dominant eigenvalue lambda 

%take the eigenvalue of maximum absolute magnitude
lambda = max( abs( diag(val) ))


%then find the stable age distribution
% by taking the corresponding right eigenvector,
% assuming that the max is in the first position
% check this with

lamsd=val(5,5)
%is this the same as above
% if so,

ssd = vec(:,5)/sum(vec(:,5))

figure(1)
bar(ssd)
xlabel('Age')
ylabel('Relative Abundance')
title('Stable Age Distribution for Caretta')
%makes a bar graph of the stable stage distribution

%to get the reproductive value vector, do the 'eig' on the transpose of the matrix
% and then follow a similar procedure

[rv_vec,rv_val] = eig(A')

%where is the dominant eigenvalue?
lamrv=rv_val(2,2)

rvs = rv_vec(:,2)/sum(rv_vec(:,2))
figure(2)
plot(rvs)
xlabel('Age')
ylabel('Relative reproductive value')
title('Reproductive Values for Caretta')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%HOW DO THESE ANALYTICAL RESULTS COMPARE WITH YESTERDAY'S
%  PROJECTION EXERCISE RESULTS????
%COMPARE THE SHAPES OF SSD, RV AND 
%COMPARE THE MAGNITUDE OF THE POPULATION GROWTH RATE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Now we will check out all the other eigenvalues and get deeper into this thing
% Make a figure of all the eigenvalues in the complex plane
figure(3)
compass(val);
xlabel('Real Axis')
ylabel('Imaginary Axis')
title('Eigenvalues of the Matrix for Caretta')

%note that the compass is normalized to the maximum lambda and then grid lines 
%  are drawn to divide the space in quarters
%  to find out the length of each arrow, take the absolute magnitude of each lambda
%investigate the eigenvalues more, take the diagonal of the d matrix

vals=diag(val)
abs(vals)
%gives you the lengths of the arrows in the Figure 4.6
%find the real parts
rs=real(vals)

%to investigate the angles in the complex plane we make use of 
% a MATLAB function that finds the arctangent with reference to
% the whole circle...  the simple "artan" command does not 
% define the angles with respect to the whole circle

anglesrad=angle(vals)

%which is in radians... to convert to degrees:

anglesdeg=anglesrad*(180/pi)

% Look at the angles in radians and deg. 
% What is the meaning of this figure?
%find the imaginary parts
is=imag(vals)

%make a figure in the complex plane the standard way, plot(x,y,'*')
%and label each eigenvalue in the complex plane
figure(4)
plot(rs, is, '*')
grid
title('Eigenvalues for Caretta, not compass')
%the real part is along the x axis and the imaginary along the y axis
%however you can also get there by angle and distance... the length of the distance is the 
%absolute magnitude, for real numbers that is just the distance from the origin along the
%x axis...but not for complex numbers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% What is the meaning of this figure?
% How is this figure different from the previous?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The next commands clear previous work
%  in order to practice the "load" command
%     and then to introduce you to a function file that finds the
%  dominant eigenvalue and corresponding left and rt eigenvectors
%  PLUS even more stuff of interest

%Practice "load" and "senstuff"
clear all
load Caretta
who
 
%use the function is called "senstuff"
%   what are the inputs?
%   what are the outputs?

%    specifying the names for the outputs
%      and providing it with an input argument

[lmat wmat vmat senmat elamat]=senstuff(A); 
echo on

vmat=vmat'
wmat

roweigen = vmat*A
coleigen = A*wmat

vmat=vmat/sum(vmat)
roweigen=roweigen/sum(roweigen)

wmat=wmat/sum(wmat)
coleigen=coleigen/sum(coleigen)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%What about complex ones?
%Find the complex eigenvalues and eigenvectors for the turtle example
%Look in "val" and "vec"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

newvec3=A*vec(:,3) % get new vector from A * vector

magvec=abs(vec(:,3)) %get absolute magnitude of original vector
magnewvec3=abs(newvec3) %get absolute magnitude of new vector

magvec=magvec./sum(magvec)  %rescale the original vector
magnewvec3=magnewvec3./sum(magnewvec3) %rescale the new vector

magval=abs(val(3,3))  %get absolute magnitude corresponding eigenvalue
angleval=angle(val(3,3)) % analyze magnitude by its angle or rotation
degrees=angleval*180/pi  % angle in degrees

figure(5)
subplot(1,2,1)
compass(vec(:,3))
title(['abs magnitude of \lambda is ',num2str(magval)])
subplot(1,2,2)
compass(newvec3)
title(['angle of rotation of \lambda is ', num2str(degrees), ' deg'])
gtext('action of the matrix on a complex eigenvector')
message=['Dont forget to manually place the text on the figure!!!!'] 
message=['Dont forget to manually place the text on the figure!!!!'] 
%don't forget to manually label the figure 

% So the matrix causes a _complex_ eigenvector to rotate 
% and changes the size of its arrows, but not their relative lengths

clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%One general thing to consider: what happens when eigenvalues of
%different kinds are raised to powers?
%Carry out exercises 3 and 4 from "creature08_eigens_day2.doc"

%Make some general rules: 
%if ?>1, 
%if ?=1, 
%if ? is positive but <1, 
%if ? is negative but >-1, 
%if ?=-1,
%if ?<-1., 
%do the powers of ? increase exponentially,
%decrease exponentially, 
%exhibit damped oscillation, 
%exhibit undamped oscillation, 
%exhibit diverging oscillation?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%What about the 
% sums of powers of conjugate pair of complex eigenvalues (Caswell, p. 78)

%consider two complex eigenvalues
l1=0.7+0.5i
l2=0.9+.6i

% find their conjugates
conl1=conj(l1)
conl2=conj(l2)

% find the sum of powers of the pair
%first set up the vector to recieve the sum
sumpairs=[ ]

%then the do loop
for t=1:40
   sumpair=l1^t + conl1^t;
   sumpairs=[sumpairs sumpair];
 end  
 
figure(6)
plot(1:40,sumpairs)
xlabel('time')
ylabel('\lambda^t + (conj\lambda)^t)')
message=['Dont forget to manually place the text on the figure!!!!'] 
gtext('\lambda = 0.7 + 0.5i')  %manually places text on the figure

%repeat for the other figure

% find the sum of powers of the pair
%first set up the vector to recieve the sum
sumpairs=[ ]

%then the do loop
for t=1:40
   sumpair=l2^t + conl2^t;
   sumpairs=[sumpairs sumpair];
 end  
 
figure(7)
plot(1:40,sumpairs)
xlabel('time')
ylabel('\lambda^t + (conj\lambda)^t)')
message=['Dont forget to manually place the text on the figure!!!!'] 
gtext('\lambda = 0.9 + 0.6i')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now look back at the projection equation for the turtles.
% Is there a sum of powers of a complex conjugate pair of eigenvalues?
% What other kinds of powers of eigenvalues are there?
% Examine the relative effects of them all over  the long term 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load Caretta
val1_ts=[];
val2_ts=[];
val34_ts=[];%complex conjugate pair sum
val5_ts=[];
for t=1:40
    val1_t=val(1,1)^t;
    val2_t=val(2,2)^t;
    val34_t=val(3,3)^t+val(4,4)^t; %complex conjugate pair sum
    val5_t=val(5,5)^t;
    val1_ts=[val1_ts val1_t];
    val2_ts=[val2_ts val2_t];
    val34_ts=[val34_ts val34_t];
    val5_ts=[val5_ts val5_t];
end

figure(8)
plot(1:40,val1_ts,'r')
hold on;
plot(1:40,val2_ts,'b')
plot(1:40,val34_ts,'c')
plot(1:40,val5_ts,'g')
forlegend=[num2str(diag(val))];
thing=[forlegend(1:3,:);
       forlegend(5,:)]
legend([thing])       
xlabel('time')
ylabel('Powers of eigenvalues')