%Calculates extinction rates of populations - using one matrix, or 
%it can draw randomly one matrix from a group in each step
%It cal also include removel of individuals from a given stage in each time step
%- e.g. simulation gardeners picking up the plants
%After multiplication of the population vector by the matrix the resulting values are replaced by a value 
%drawn from Poisson distribution
%also the number of removed plants is not fixed but replaced by a value drawn from a Poisson distribution 
%with a given mean
function gardeners (mat,k,tlimit,m,groups,vector,rep,zah,nozah);

% zahradkari('mat.txt',3,10,2,2,'vec.txt',8,'zah.txt',3);
% cd c:\docasne\skripty-td, pridat poisson k nasobeni i zahradkarum

% k - number of stages
% tlimit - number of years
% m - number of matrices per group
% groups - no. of groups of matrices, vector=vector of initial number of individuals per stage
% rep - number of replicates
% zah - row vector containing number of individuals sampled per per year from a given stage, threre can be more rows - more versions of sampling 
% nozah - number of rows in file zah 

pois=1;
fid=fopen(mat,'r'); %read the matrices
b=fscanf(fid,'%f',[k,k*m*groups]);
l=0;
for i=1:m*groups
    for j=1:k
      for p=1:k
   	    aa (j,p,i)= b(j,p+k*(i-1));
       end
   end
aa(:,:,i)=transpose(aa(:,:,i));   
end
fclose(fid);

fid=fopen(vector,'r');
vecc=fscanf(fid,'%f',[k*groups,1]);
fclose(fid);
fid=fopen(zah,'r');
zah=fscanf(fid,'%f',[k,nozah]);
fclose(fid);
for ia=1:groups % groups of matrices
   for ja=1:k
      vec1(ja,1)=vecc(ja+k*(ia-1));
   end   
   for ja=1:m
      a(:,:,ja)=aa(:,:,(ja+(ia-1)*m));
   end
   for noz=1:nozah % gardeners run
       for r=1:rep 
           vec=zeros(k,tlimit+1);
           vec(:,1)=vec1;
           % generate sequence of environments
           env=ceil(rand(1,tlimit)*m);
           %iteration
           for t=1:tlimit % time
               c=a(:,:,env(t));
               vec(:,t+1) = c*vec(:,t);
               if pois==1
               for tt=1:k
                  vec(tt,t+1)=poisr(vec(tt,t+1));
                  zzzah(tt,noz)=poisr(zah(tt,noz));
               end
               else
                  zzzah=zah;
               end   
               vec(:,t+1)=vec(:,t+1)-zzzah(:,noz);
               for u=1:k
                  if vec(u,t+1)<0
                     vec(u,t+1)=0;
                  end
               end   
               suma(t+1)=sum(vec(:,t+1));
               fid5=fopen('popsize.out','a');
               if ia==1 & noz==1 & r==1 & t==1
                  fprintf(fid5,'time\tgroup\trep\tzah\tpopsize\n');
               end   
               fprintf(fid5,'%10.7f\t%10.7f\t%10.7f\t%10.7f\t%10.7f\t',t+1,ia,r,noz,suma(t+1));
               for u=1:k
                  fprintf(fid5,'%10.7f\t',vec(u,t+1));
                  if u==k
                  fprintf(fid5,'\n');
                  end
               end
               fclose(fid5);
           end
       end
    end
 end    

 function y = poisr(lambda) %,n)
% Uses naive inversion method.
if lambda == 0
   for i = 1 %:n,
      y(i) = 0;
   end;
   return;
end

if nargin==1, n=1; end;
low = max( 0 , floor( lambda - 8*sqrt(lambda) ));
hi = ceil( lambda + 8*sqrt(lambda) + 4/lambda );
x = (low:hi)';
p = cumsum( exp( -lambda + x.*log(lambda) - gammaln(x+1) ));
u = rand(n,1);
y = zeros(n,1);
for i = 1%:n,
   k = find( u(i) < p );
   y(i) = x( k(1) );
end