% Simple methods for calculating age-based life history parameters for 
% stage-structured populations
% Cochran et Ellner (1992), Ecol. Monogr. 62(3): 345-364.
function cochran(Amat,Bmat,Pmat,Dmat, stg, time,repro)
% Amat = original matrix of the species 
% Bmat = matrix with transitions of reproduction
% Pmat = matrix with transitions of growth and survival
% stg = number of stages
% time = number of years until which survivorship probability and stable age 
% distribution should be calculated
% Dmat = set of matrices Pmat, one bellow each other,  with all transition in 1st, 2nd.... 
% column equal to zero
% i.e. - all individuals entering stage 1,2... will die in the next year
% repro = matrix with columns meaning reproduction contaning only zeros - reproductive individuals
% do not survive

%reads Amat
fid1=fopen(Amat,'r');
A=fscanf(fid1,'%f',[stg,stg]);
A=transpose(A);
fclose(fid1);

%reads Bmat
fid1=fopen(Bmat,'r');
B=fscanf(fid1,'%f',[stg,stg]);
B=transpose(B);
fclose(fid1);

%reads Pmat
fid1=fopen(Pmat,'r');
P=fscanf(fid1,'%f',[stg,stg]);
P=transpose(P);
fclose(fid1);

%inictializes survivorship function
l=zeros(time,stg);

%reads Dmat
fid1=fopen(Dmat,'r');
DD=fscanf(fid1,'%f',[stg,stg*stg]);
fclose(fid1);
D=zeros(stg,stg,stg);
for i=1:stg
   for j=1:stg
      for k=1:stg
         D(j,k,i)=DD(j,k+stg*(i-1));
      end
   end
end   
for i=1:stg
   D(:,:,i)=transpose(D(:,:,i));
end   

%reads repro
fid1=fopen(repro,'r');
Q=fscanf(fid1,'%f',[stg,stg]);
fclose(fid1);
Q=transpose(Q);

% survivorship function
for j = 1:stg
   for t = 2:time
	    for i = 1:stg     
           l(t,j)= l(t,j)+P(i,j)^(t-1);
       end
   end
end 
%print results of the survivorship function
fid=fopen('surfunc.out','w');
fprintf(fid,'stage\t expected time\t probability\n');
fclose(fid);
for j=1:stg
	 for t=2:time   
		  fid=fopen('surfunc.out','a');
		  fprintf(fid,'%10.7f	%10.7f	%10.7f\n',j,t,l(t,j));
        fclose(fid);
    end
 end    
 
 % remaining life span
 % generate identity matrix
 I=zeros(stg,stg);
 for i=1:stg
    I(i,i)=I(i,i)+1;
 end
 % end of generation of identity matrix
 L=zeros(stg,1);
 minus=(I-P)^-1;
 for j=1:stg
    for i=1:stg
        L(j)=L(j)+minus(i,j);
    end
 end
% print results of the remaining life span
fid=fopen('remain.out','w');
fprintf(fid,'stage\t remaining life span\n');
fclose(fid);
for j=1:stg
    fid=fopen('remain.out','a');
	 fprintf(fid,'%10.7f\t %10.7f\n',j,L(j));
    fclose(fid);
end    

% conditional total life span of newborn of type i, 
% conditional on passing through stage j at least once
for i=1:stg
    D(:,:,i);   
    minus(:,:,i)=(I-D(:,:,i))^-1;
    minus2(:,:,i)=(I-D(:,:,i))^-2;
    for j=1:stg
       Etau(i,j)=minus2(i,j)/minus(i,j);
    end
 end
 for i=1:stg
	  for j=1:stg 	    
         LAMBDA(i,j)=Etau(i,j)+L(i);
      end
 end     
% print results of the conditional life span
fid=fopen('conditional.out','w');
fprintf(fid,'newborn type-j\t stage-pass-i\t life span\n');
fclose(fid);
for i=1:stg
	 for j=1:stg   
    	  fid=fopen('conditional.out','a');
	     fprintf(fid,'%10.7f\t %10.7f\t %10.7f\n',j,i,LAMBDA(i,j));
        fclose(fid);
    end
end     

%age at first reproduction
minus=(I-Q)^-1;
minus2=(I-Q)^-2;
Ealpha=zeros(stg,1);
for j=1:stg
    for i=3 %code of reproductive stage  
        Ealpha(j)=Ealpha(j)+minus2(i,j)/minus(i,j);
    end   
end
% print results for age at first reproduction
fid=fopen('first repro.out','w');
fprintf(fid,'stage\t first repro\n');
fclose(fid);
for j=1:stg
    fid=fopen('first repro.out','a');
	 fprintf(fid,'%10.7f\t %10.7f\n',j,Ealpha(j));
    fclose(fid);
end    

%population averages of life history traits
[W,d] = eig(A);
lambda = diag(d);
imax = find(lambda==max(lambda));
V=conj(inv(W));
w = W(:,imax);
v=real(V(imax,:));
lambda1 = lambda(imax);
Bw=B*w;
for j=1:stg
    b(j)=Bw(j)/sum(Bw);
end  
 
%Population dynamics 
%mean age in stage classes
one=(I-((lambda1^-1)*P))^-1;
two=(I-((lambda1^-1)*P))^-2;

for i=1:stg
    up=0;
    down=0;   
    for j=1:stg    
        up=up+two(i,j)*b(j);
        down=down+one(i,j)*b(j);
    end
    y(i)=up/down; 
end    

% print results for mean age in stage classes 
fid=fopen('mean age stage.out','w');
fprintf(fid,'stage\t mean age\n');
fclose(fid);
for i=1:stg
    fid=fopen('mean age stage.out','a');
	 fprintf(fid,'%10.7f\t %10.7f\n',i,y(i));
    fclose(fid);
end    

%age composition of stage classes
% be careful - it does not sum to one!!!!!  
for i=1:stg
    for t=1:time
        ca=P^(t-1);
        up=0;
        down=0;
        for j=1:stg
            up=up+ca(i,j)*b(j);
            down=down+one(i,j)*b(j);
        end    
        p(i,t)=((lambda1^(-t-1))*up)/down;
    end
end
  
% print results for age composition of stage classes
fid=fopen('age composition by stage.out','w');
fprintf(fid,'stage\t age \t fraction\n');
fclose(fid);
for i=1:stg
    for t=1:time
        fid=fopen('age composition by stage.out','a');
	     fprintf(fid,'%10.7f\t %10.7f\t %10.7f\n',i,t,p(i,t));
        fclose(fid);
    end
end     

%stable age distribution
for t=1:time
    up=0;   
    down=0;
    for i=1:stg
        up=up+p(i,t)*w(i);
        down=down+w(i);
     end
     omega(t)=up/down;
end     

% print results for stable age distribution
% be careful - does not sum to one
fid=fopen('stable age distribution.out','w');
fprintf(fid,'age \t fraction\n');
fclose(fid);
for t=1:time
    fid=fopen('stable age distribution.out','a');
	 fprintf(fid,'%10.7f\t %10.7f\n',t,omega(t));
    fclose(fid);
end     

% Maternity function  
gama=zeros(stg);
for i=1:stg
    for j=1:stg
        gama(i)=gama(i)+B(j,i)*(v(j)/v(1));
    end   
end     
for t=1:time
     pt=P^(t-1);
     for j=1:stg
         up=0;
         down=0;
         for i=1:stg        
             up=up+pt(i,j)*gama(i);
             down=down+pt(i,j);
         end   
         f(t,j)= up/down;
     end
end
   
% print results for maternity function
fid=fopen('maternity function.out','w');
fprintf(fid,'stage \t age\t offspring production measured in sgt1 values \n');
fclose(fid);
for j=1:stg
    for t=1:time
        fid=fopen('maternity function.out','a');
	     fprintf(fid,'%10.7f\t %10.7f\t %10.7f\n',j,t,f(t,j));
        fclose(fid);
    end
end     

% Net reproductive rate
Ro=zeros(j);
for j=1:stg
    for t=1:time
        Ro(j)=Ro(j)+l(t,j)*f(t,j);
    end
end    

% print results for net reproductive rate
fid=fopen('net reproductive rate.out','w');
fprintf(fid,'stage \t net reproductive rate\n');
fclose(fid);
for j=1:stg
    fid=fopen('net reproductive rate.out','a');
	 fprintf(fid,'%10.7f\t %10.7f\n',j,Ro(j));
    fclose(fid);
end     

%generation time
ic2=(I-P)^(-2);
for j=1:stg
    up=0;
    for i=1:stg   
        up=up+ic2(i,j)*gama(i);
    end
    mi(j)=up/Ro(j);
end    

% print results for generation time
fid=fopen('generation time.out','w');
fprintf(fid,'stage \t generation time\n');
fclose(fid);
for j=1:stg
    fid=fopen('generation time.out','a');
	 fprintf(fid,'%10.7f\t %10.7f\n',j,mi(j));
    fclose(fid);
end