% PROGRAM harvestopt.m
% This MATLAB code simulates a simple population growth situation 
% in which a population increases by one individual each time period 
% with probability p and decreases by one individual each time 
% period with probability 1-p. The population is harvested when it
% reaches level b, and a fraction f of the population is harvested
% at that time, so that the population level then drops to the 
% greatest integer less than b*(1-f). This code allows the user to
% investigate how large b and f can be made, for particular values of p, 
% while maintaining the population above a certain minimum acceptable
% level with high probability (this is the minimum viable population
% size).  
%  Thus, this is the level-crossing
% calculation in which one wants to choose the level b at which to 
% allow harvesting to occur so as to maximize the harvest while 
% ensuring with high probability that the population level does not
% fall too low. The population is assumed to start with 100 individuals.
% The minimum viable population size is assumed to be .1 of this
% that is 10 individuals. 
%     The user is prompted to supply
% the probability p, the number n of simulations for each b value 
% used, the farction f of the population harvested and the range of b-values
% to use (the code partitions this range into mm equal pieces). 
% At the end, output will be a vector, x, giving the number of harvests
% and a vector y giving the total amount of harvest 
% for each simulation. The total simulation time is 300 time units.  
% The mean and standard deviation of each of these vectors is saved for
% each b value in the vectors mtim,stim for x, and mtotharv, stotharv for y. 
%     Graphs of these means and std are plotted as a function of the
% b values.
%
p=input('p = probability of going up one unit:     ')
f=input('f = fraction of population harvested when level b is reached:     ')
n=input('n= number of simulations for each b:'   )
blow=input('blow=lower value for range of b at which to harvest:   ')
bhi=input('bhi=upper value for range of b at which to harvest:   ')
x=[ ];
y=[ ];
c=[ ];
mtotharv=[ ];
stotharv=[ ];
mtim=[ ];
stim=[ ];
popinit=100;
minviable=.1*popinit;
% mm is the number of cases for b to run
mm=10;
bval(1)=blow;
for i=2:(mm+1)
   bval(i)=blow+(i-1)*(bhi-blow)/mm;
   end;
% iterate over the mm+1 b values
bval
for ii=1:(mm+1)
bb=bval(ii);
% for each b value, make n simulations
% here the current state is xn, total number of harvests is
% ttot and total harvest is totharv
for j=1:n
ttot=0;
totharv=0;
itime=0;
xn=popinit;
    while (itime < 300);
           if (xn < minviable)
             ttot=0;totharv=0; 
             break;
           else
           end;
       if (xn < bb)
         a=rand(1,1);
         if a <= p
           xn=xn+1;
         else 
           xn=xn-1;
         end;
         itime=itime+1;
       else
          ttot=ttot+1;
          totharv=totharv+bb*f;
          xn=bb*(1-f);
        end;
     end;
% The vector x contains the total number of harvests
% the vector y contains the total harvest for each simulation 
 x=[x ttot];
 y=[y totharv];
end;
mtim=[mtim mean(x)];
stim=[stim std(x)];
mtotharv=[mtotharv mean(y)];
stotharv=[stotharv std(y)];
x=[ ];
y=[ ];
end;
plot(bval,mtim,'o',bval,stim,'+')
title('mean number of harvests for various b values as o, std as +'),pause
plot(bval,mtotharv,'o',bval,stotharv,'+')
title('mean total harvest for various b values as o, std as +'),pause
disp('Values of harvest level b used')
bval,pause
disp('Mean number of harvests for each harvest level')
mtim,pause
disp('Mean total harvest for each harvest level')
mtotharv,pause