clear; close all

% Script for running BayesForaG and related functions. When this is run a
% Bayesian policy is generated. This is set up here so that it finds the
% long-term intake rate maximizing policy, but the step where the Bayesian
% strategy is found could be set up to maximize any other fitness currency.
% It uses a negative binomial prey distribution, but any other prey
% could be used.
% 
% Three subplots are created in one figure, to show some of the inputs and
% outputs generated.
% 
% See also BayesForaG, Sim_nt, NegBin

% Ola Olsson, 2009-05-31

%% Initializing parameters 
nmax=60;   % Max number of prey to consider in the prey distr. used by NegBin
C=6;       % Initial guess of long-term intake rate
A=1;       % Searching efficiency
T=2;       % Maximum search time
D=.05;     % Time step
k=100;     % Number of patches of each quality that are generated in the simulation
tau=.01;   % Travel time

%% Initializing prey distribution
% negative binomial prey distribution generated here, but any discrete
% distribution works
[nmax, pN]=NegBin(6, 1, nmax);  
figure(1); subplot(3, 1, 1)
bar(0:nmax, pN)
ax=axis;
axis([-.5 nmax+.5 ax(3:4)])
xlabel('Initial prey density')
ylabel('Frequency')


%% Initializing cPn
%  See Olsson & Brown (in press, Oikos) and Sim_nt for explanation.
tic
Pn=zeros(nmax+1, nmax+1);
for Ni=0:nmax
    n=0:Ni;
    Pn(n+1, Ni+1)=PnT((0:Ni)', D, Ni, A);
end
cPn=cumsum(Pn);
toc

%% Finding the best Bayesian strategy
% The strategy found here is simple rate-maximization as we find
% C*=nbar/(tbar+tau), i.e. the aimed for C and the actual long-term rate
% are the same.

d=1; i=0;
while d>.01
    i=i+1;
    tic
    % This step finds the stopping points for the prey distribution and C
    [SP, X, G, S]=BayesForaG(pN, C, A, T, D);         
    
    % This step simulates the prey capture in k patches of each N
    [nbar, tbar, n, t, GUD, N, pNx]=Sim_nt(pN, X, D, cPn, k);
    elt=toc;
    C0=C;
    C=nbar./(tbar+tau);   % This is the long-term intake rate using the present strategy
    d=abs(C0-C);          % Difference between previous and this strategy
    disp([i C d elt])
end

%% Plotting the best strategy
subplot(3, 1, 2)
nx=[0 nmax]; tx=[0 T];
imagesc(tx, nx, X)
set(gca, 'YDir', 'Normal')
xlabel('Search time, t')
ylabel('Prey caught, n')

%% Plotting GUD by search time
% As k patches of each quality are simulated simple mean values cannot be
% used, but weighted means need to be calculated as below
subplot(3, 1, 3); hold on
tu=unique(t);
for ti=tu'
    ix=find(t==ti);
    mGUD=sum(pNx(ix).*GUD(ix))./sum(pNx(ix));
    plot(ti, mGUD, 'o')
end
xlabel('Search time')
ylabel('Giving-up density, GUD')