function [SP, X, G, S]=BayesForaG(pN, C, A, T, D)

% [SP, X, G, S]=BayesForaG(pN, C, A, T, D)
%
%   A model for Bayesian foraging as described in Olsson & Brown (2009).
%   The G in the function name stands for 'general' or for 'Green', as in
%   Dick (Richard F.) Green, who has inspired this work for many years.
% 
%   You are free to download and use this code in your own work, on two
%   simple conditions. First, please, cite the relevant publictations when
%   using this for your own publications. Papers to cite are listed below.
%   Second, do send me an email to let me know that you are using the
%   stuff: ola.olsson@zooekol.lu.se . With some luck I may also be able to
%   provide some help, if needed. Also, please report any bugs or mistakes
%   you find.
%   
%   The model handles only a forager using random search, but any prey
%   distribution.
%
%   The function returns a vector of optimal stopping points, SP, a policy 
%   matrix, X, a matrix of expected gains, G, and a matrix of expected
%   times, S, to remain in the patch. SP or X is normally the usueful 
%   output.
% 
%   The function requires five inputs. pN is a vector of the frequency 
%   ditribution of prey densities in patches. It is assumed that the first
%   element in pN corresponds to 0 prey items, the second to 1 item, and so
%   on. The length of the vector will be checked to calculate Nmax. There
%   is no limitation on Nmax, but the computing times required for large
%   Nmax (say Nmax>50) increases rapidly. C is the potential inatake rate
%   at which the patch should be abandoned. A is the searching efficiency.
%   The model has not been rigorously tested for values other than A=1, but
%   should work for any value. T is the maximum allowed search time. Often
%   1.5>T>4 is sufficient when A=1. D is the time step in the model.
%
%   The idea to abandon patches based on the ratio of G/S was first
%   suggested by Green in his early papers (e.g. Green, 1980) and is
%   related to McNamara's (1982) potential function. See Green 2006, or
%   McNamara et al. 2006 for recent reviews of this.
%
%   Green, R. F. 1980. Bayesian birds: a simple example of Oaten's
%       stochastic model of optimal foraging. - Theor. Popul. Biol. 18: 244-256
% 
%   Green, R. F. 2006. A simpler, more general method of finding the
%      optimal foraging strategy for Bayesian birds. - Oikos 112: 274-284.
% 
%   McNamara, J. M. 1982. Optimal patch use in a stochastic environment. -
%       Theor. Popul. Biol. 21: 269-288.
% 
%   McNamara, J. M., Green, R. F. and Olsson, O. 2006. Bayes' theorem and
%      its applications in animal behaviour. - Oikos 112: 243-251.
% 
%   Olsson, O & Holmgren NMA. 1998. The survival-rate-maximizing policy for
%       Bayesian foragers: wait for good news. Behav.Ecol. 9: 345-353.
% 
%   Olsson, O & Brown, J. S. In press. Smart, smarter, smartest: Foraging
%      information states and coexistence. Oikos.
% 
% See also XoverN, PnT, Posterior
 
% OO - 2008-08-13

warning off MATLAB:divideByZero % The calculation of r towards the bottom often involves 0 in the denominator, which is no error

% Adjustments that should solve some common issues with prey distributions
ix=find(pN); pN=pN(1:max(ix));  %#ok<MXFND> 
Nmax=length(pN)-1;  
pN=pN./sum(pN);

tx=0:D:T;
tix=1:length(tx);
nx=(0:Nmax)';
G=zeros(length(nx), length(tx)); S=G; 
komb=XoverN(Nmax);        % Binomial coeffs up to Nmax
Pk=zeros(Nmax+1, Nmax+1);
for x=0:Nmax  % prey left in the patch
    k=0:x;    % prey to catch in next time step
    Pk(x+1, k+1)=PnT(k, D, x, A);  % probability to catch k prey in next time step given x prey left. O&H 2000
end

for ti=fliplr(tix(1:end-1))
    t1=tx(ti);
    Rx=Posterior(Nmax, t1, pN, A, komb);   % Posterior probability see the func for explanation
    for n=0:Nmax % Number of prey taken so far    
        sk=zeros(Nmax-n, Nmax-n);
        for N=n:Nmax   % N that could have been in patch initially, given that n are taken now
            x=N-n;
            for k=0:x
                sk(k+1, x+1)=Pk(x+1, k+1).*Rx(n+1, x+1); % see Sk below
            end
        end
        Sk=sum(sk, 2);  % probability of catching k items in next time step. Eq. 3 in O&H 1998, Eq. 7 in O&B in press
        kx=(0:(length(Sk)-1))'; n2=n+kx;
        G0=Sk.*(kx+G(n2+1, ti+1));  % see G below
        G(n+1, ti)=sum(G0);         % expected number of prey to catch from now until the patch is left. Eq. 9 in O&B 2009, see also Green 2006
        
        S0=Sk.*(D+S(n2+1, ti+1));   % see S below
        S(n+1, ti)=sum(S0);         % expected remaining time to spend in patch from now on. Eq. 10 in O & B 2009. See also Green 2006
    end
    r=G(:, ti)./S(:, ti);           % potential intake rate; the rule to base patch departure on; the holy grail; see O&H 1998, O&B 2006 or Green 2006 for details
    ix=find(r<C); G(ix, ti)=0; S(ix, ti)=0;    % patches with potential intake rate less than aimed-for are abandoned
end
X=G>0;  % The policy as 0 and 1
SP=sum(X, 2).*D;  % The stopping points. This may be messed up for strange distributions. X is always correct.
warning on MATLAB:divideByZero