% This is a script that can boostrap various parts of the matrix - the raw data, but also diferent tyes of fixed data
% provided in basmat and pdpger files
% it works with only one transition interval
function boot4(file,noboot,stg,rows,basmat,pdpger,confidence,plot,noind,Gnoind,reprodukce)
% file - file with data containing 2 + reprodukce columns - stage year 1, stage year 2, seed production into stage 1, into stage 2...
% Each row of data has information on one individual.
% noboot - number of bootstrap values - 1, give 0 if bootstrap not required
% stg - number of stages
% rows - number of individuals = number of rows  
% Basmat contains all data not based on file it is added to the resulting matrix 
% basmat is multiplied by nonzero transitions in pdpger
% pdpger - matrix with probability of germination of seed for each stage, contains zero for transitions that are not germination, 
% nonzero elements of pdpger are used to multiply the resulting matrix
% confidence is alpha for ci
% plot the distribution of the resulting lambda
% noind - vector of number of individuals for bootstrap of basmat-it is
% possible to bootsrap only values below 1
% Gnoind - vector of number of individuals for bootstrap of pdpger
% NOTE that transions in basmat larger than 1 are not nebootstraped
% reprodukce - number of columns with reproduction in file - used when reproduction involves not only individuals into stage 1 but also to other stages

% reads basic matrix:
notrans=1; % do not change!!!
fid1=fopen(basmat,'r'); % reads basmat
b=fscanf(fid1,'%f',[stg,stg*notrans]);
b=transpose(b);
fclose(fid1);
cols=2+reprodukce;
if strcmp(file,'none') ~= 1 % if file exists
%reads matrix with data
	fid=fopen(file,'r');
	a=fscanf(fid,'%f',[cols,rows]);
	a=transpose(a);
	fclose(fid);
%read germination probability matrix
	fid2=fopen(pdpger,'r');
   c=fscanf(fid2,'%f',[stg,stg]); 
   c=transpose(c);
	fclose(fid2);
end
% reads noind in basmat
if strcmp(noind,'none') ~= 1
	fid2=fopen(noind,'r');
	noind=fscanf(fid2,'%f',[1,stg]);
	noind=transpose(noind);
   fclose(fid2);  
end

% reads Gnoind = noind in pdpger
if strcmp(Gnoind,'none') ~= 1
	fid2=fopen(Gnoind,'r');
	Gnoind=fscanf(fid2,'%f',[1,stg]);
	Gnoind=transpose(Gnoind);
   fclose(fid2);  
end

lam = zeros(size(noboot+1));%initiates zero matrices of the required size
ela = zeros (stg,stg, noboot+1); %size of the matrix,size of the matrix,no of replicates
sen = zeros (stg,stg, noboot+1); %size of the matrix,size of the matrix,no of replicates
test2 = zeros(rows,cols);%initiates zero matrices for bootstrapped values
%function bootstrap (noboot, stg,notrans,a,rows, notrans) 
%start of the bootstrapping part
%strcmp(noind,'none')
for ee = 1:noboot+1 %number of replicates
   
   for z=1:stg
   if sum(sum(b))>0 & (noind(z))~=0 & notrans==1 & noboot>0 & max(b(:,z))<=1% if basmat contains data, if noind is nonzero
       b;
       stg;
       noind;
       reprodukce;
       z;
       bb=bootbas(b,stg,noind,reprodukce,z);
       kuk=bb;
         bas=bootfile(noboot,bb,noind(z));
         gg=1;
         hh=ones(1,stg);
         b0=makemat(bas,notrans,stg,noind(z),hh,gg,reprodukce);
         b1(:,z)=b0(:,z);
         gg=0;
      else
         b1(:,z)=b(:,z);  %bootbas (b,stg,noind);    
    end
    end
    
    % bootstrap of pdpger
    for z=1:stg
    if sum(sum(c))>0 & Gnoind(z)>0 & notrans==1 & noboot>0 % if pdpger contains data, if noind is nonzero
	      cc=bootbas (c,stg,Gnoind,reprodukce,z);
         ccc=bootfile(noboot,cc,Gnoind(z));
         gg=1;
         hh=ones(1,stg);
         c0=makemat(ccc,notrans,stg,Gnoind(z),hh,gg,reprodukce);
         c1(:,z)=c0(:,z);
         gg=0;
    else c1(:,z)=c(:,z);
    end 
    end

	 if strcmp(noind,'none') == 0 % if bootstraping basmat
       test2=bootfile(noboot,a,rows);
       for y=1:stg
	   	  for z=1:stg	            
               if b(y,z)==0 | b(y,z)>1 | noind(z)==0
                  b1(y,z)=b(y,z);
               end
            end
        end    
        if noboot == 0 
      	  test2 = a;
        end % end of the bootstraping part
        gg=0;
        tran = makemat(test2,notrans,stg,rows,c1,gg,reprodukce);
        mat = addbasmat(notrans,stg,b1,tran,c1,c);
    end
	 if strcmp(file,'none') == 1
		 mat=b1
	 end
    
    for i=1:notrans   
      for y=1:stg
             for yy=1:stg
   fid14=fopen('matice.out','a');
   fprintf(fid14,'%10.7f\t%',mat(y,yy,i));
   if yy==stg
      fprintf(fid14,'%s\t%s\t%10.7f\t%10.7f\n',file,basmat,ee,y);
   end   
   fclose(fid14);
   end
   end
      
      [W,d] = eig(mat(:,:,i));
		lambda = diag(d);
		imax = find(lambda==max(lambda));
		V=conj(inv(W));
		w = W(:,imax);
		v=real(V(imax,:));
		lambda1 = lambda(imax);
		sen = v'*w';
		ela = sen.*mat(:,:,i)/max(eig(mat(:,:,i)));
		one = size (lambda1);
		senM(:,:,i,ee) = sen;
		elaM(:,:,i,ee) = ela;
      
     if stg==3
         for y=1:3
   fid14=fopen('ela.out','a');
	fprintf(fid14,'%10.7f\t%10.7f\t%10.7f\t%s\t%s\t%10.7f\n',elaM(y,1,i,ee),elaM(y,2,i,ee),elaM(y,3,i,ee),file,basmat,ee);
   fclose(fid14);
         end
      end
      if one (1) >1
  			lambda1 = lambda1(1);
		end
		lam(ee,i) = lambda1;
	end
end
if noboot > 99      
   if plot==1
      figure;
      hist(lam(:,i));
      TITLE (sprintf ('transition number: %1d', notrans));
   end
	ci = zeros(notrans,2);
	ciS = zeros (stg*notrans,stg*2);
	ciE = zeros (stg*notrans,stg*2);
	k=floor((noboot+1)*confidence./2);
	lowind= k;
	highind=noboot - k +1;
  	for i = 1:notrans
     for j = 1:stg
         for k = 1:stg
            o=0;
           	for ee = 1:noboot
               o = o+1;
               senMM (o)= senM (j,k,i,ee);
               elaMM (o)= elaM (j,k,i,ee);
            end 
            sortedS =sort(senMM);
            ciS(j+stg*(i-1),(2*k)-1)=sortedS(lowind);
            ciS(j+stg*(i-1),(2*k))= sortedS(highind);
            sortedE = sort(elaMM);
     			ciE(j+stg*(i-1),(2*k)-1)=sortedE(lowind);
            ciE(j+stg*(i-1),(2*k))= sortedE(highind);
            bootstrap_mean = mean(lam(:,i));
   	      sorted=sort(lam(:,i));  
            ci(i,:)=[sorted(lowind),sorted(highind)];
         end
     end
  end 
	fid4=fopen('CI','a');
	fprintf(fid4,'%10.7f\t%10.7f\t%s\t%s\n',ci,file,basmat);
	fclose(fid4);
	fid6=fopen('ciS','a');
	ciS = transpose(ciS); 
	fprintf(fid6,'%10.7f\t%10.7f\n',ciS);
	fclose(fid6);
	fid7=fopen('ciE','a');
	ciE = transpose(ciE); 
	fprintf(fid7,'%10.7f\t%10.7f\n',ciE);
   fclose(fid7);
end
if noboot>0
   fid8=fopen('matice.out','a');
   for i=1:stg
      for j=1:stg
         fprintf(fid8,'%10.7f\t',mat(i,j));
         if j==stg
            fprintf(fid8,'\n');
         end
      end
   end   
   fclose(fid8);
     
end
if noboot == 0
    fid8=fopen('matice.out','a');
   for i=1:stg
      for j=1:stg
         fprintf(fid8,'%10.7f\t',mat(i,j));
         if j==stg
            fprintf(fid8,'\n');
         end
      end
   end   
   fclose(fid8);
        
   fid5=fopen('lambda','a');
	fprintf(fid5,'%10.7f\t',lam);
	fclose(fid5);
   fid6=fopen('sen','a');
	sen = transpose(sen); 
	fprintf(fid6,'%10.7f\t',sen);
	fclose(fid6);
	fid7=fopen('ela','a');
	ela = transpose(ela); 
	fprintf(fid7,'%10.7f\t',ela);
   fclose(fid7);
   fid8=fopen('mat','a');
	mat = transpose(mat); 
	fprintf(fid8,'%10.7f\t',mat);
	fclose(fid8);
end

function [tran]=makemat(test2,notrans,stg,rows,c1,gg,reprodukce)   
rep = zeros (stg,stg); % vector of no of individuals per class
no = zeros (notrans,stg);   
t = zeros(stg*notrans,stg); %matrix with no ind in each transition
for i = 1:rows %no of rows in the data file
    for j = 1:notrans
        no(j,test2(i,j)) = no(j,test2(i,j)) + 1; % counts number of individuals per stage in the 1st year
        for juj=1:reprodukce % moznost mit vic sloupcu s produkce
           if test2(i,j)~=0
              test2;
           rep(juj,test2(i,1)) = rep(juj,test2(i,1))+ test2(i,2+juj);%counts no of fruits    
        end   
        end
        if test2(i,1)~=0 & test2(i,j+1)~=0 % base to create the matrix
         	t(test2(i,2),test2(i,1)) = t(test2(i,2),test2(i,1)) + 1; % counts number of individuals per stage in the 1st year
		  end
     end
 end
 %calculates transition probabilities 
 for i = 1:stg %no of stages
     for j = 1:stg %no of stages
        for k = 1:notrans
            if no (k,j)>0
              	tran (i,j) = t(i,j) ./no(1,j);%calculates transition probabilities
            else 
           		tran (i+stg*(k-1),j) = 0;   
           	end   
        end
	  end
  end
  %adds seed production x germination
  if gg==0
     %for o = 1:notrans
     for mm=1:stg
        for m = 1:stg %no of stages
           if no(1,m)>0
              tran (mm,m)= tran (mm,m)+ (rep(mm,m)/no(1,m))*c1(mm,m);
           end   
     	   end
     end
  end      
  return;
  function [mat] = addbasmat(notrans,stg,b1,tran,c1,c)
  %adds data to basmat matrix with transition probabilities and saves matrix in a speciefic format
    for i = 1:stg %no of stages
      for j = 1:stg %no of stages
         for k = 1:notrans
            if b1(i+stg*(k-1),j)~= 0
               if c (i,j)>0
                  b1(i,j)=b1(i,j)*c1 (i,j); % multiply basmat with pdpger
               end   
                  tran (i+stg*(k-1),j)= tran (i+stg*(k-1),j)+ b1 (i+stg*(k-1),j); % Basmat se pricita!!!!
            end   
           	mat(i,j,k)=tran (i+stg*(k-1),j);
         end
     end
  end
  return;      

function [bb]=bootbas(b,stg,noind,reprodukce,z) % bootstrap of basmat, z=column
zz=0;
xx=0;
%for z=1:stg % for each stage
   xx=zz;
    for y=1:stg
      if b(y,z)~=0 & b(y,z)<=1 & noind(z)~=0 % for transitions below 1, for more then 0 noind
         for x=1:ceil(noind(z)*b(y,z))
             xx=xx+1;
             bb(xx,1)= z;
             bb(xx,2)=y;
             for t=1:reprodukce
                 bb(xx,2+t)=0;
              end    
         end     
         for hh=xx+1:noind(z)+zz; % assumes that nonzero transition cannot be below zero transition in the same column
             bb(hh,1)= z;
             bb(hh,2)=0;
              for t=1:reprodukce
                 bb(xx,2+t)=0;
              end    
         end
      end      
   %end
   zz=zz+noind(z);
end
return;

function [fileout] = bootfile(noboot,filein,size) 
if noboot > 0      
  	ii = 0;
	while ii<size %no of lines to be sampled
      ii=ii+1; 
      %size
      x=round(rand*size+0.5);
      x;
      filein;
      fileout(ii,:)=filein(x,:);
   end
else fileout = filein; 
end
return;