function X = MEBPincr(X, n, offspringftn, weightsftn)

% function X = MEBPincr(X, n, offspringftn, weightsftn)
%
% increments the state of an MEBP process
%
% X in is the current state of the process at step k (levels 0 to nmax)
% n is the level being incremented
% offspringftn is a function handle for a function that samples from 
%   the offspring distribution p(n)
% weightsftn is a function handle for a function that samples from 
%   the weights distribution, given the number of offspring
% X out is the state of the process at step k+1
%
% item n+1 of X in is (kappa(0,n,k), S(0,n,k), Z(0,n+1,k), alpha(0,n,k), P(0,n+1,k))
% item n+1 of X out is (kappa(0,n,k+1), S(0,n,k+1), Z(0,n+1,k+1), alpha(0,n,k+1), P(0,n+1,k+1))
%
% types are 0+ 0- 1+ 1- -1+ -1-, labelled 1 2 3 4 5 6
%
% Owen Jones 6 July 2003, 17 January 2007

if (n > 0) & ~(X{n}{2} == X{n}{3})
    error('X is not at end of level n crossing')
end

if X{n+1}{2} == X{n+1}{3}
    X = MEBPincr(X, n+1, offspringftn, weightsftn);
    X{n+1}{2} = 1;
    X{n+1}{3} = feval(offspringftn, []);
    X{n+1}{5} = feval(weightsftn, X{n+1}{3});
else
    X{n+1}{2} = X{n+1}{2} + 1;
end

X{n+1}{1} = X{n+1}{1} + 1;

nmax = length(X) - 1;
if n == nmax
    if X{n+1}{2} == X{n+1}{3}
        if X{n+1}{4} == 1, X{n+1}{4} = 3; else X{n+1}{4} = 6; end
    elseif mod(X{n+1}{2},2) == 1
        if rand < 0.5, X{n+1}{4} = 1; else X{n+1}{4} = 2; end
    else
        if X{n+1}{4} == 1, X{n+1}{4} = 4; else X{n+1}{4} = 5; end
    end
else
    if X{n+1}{2} == X{n+1}{3}
        if (X{n+2}{4} == 1) | (X{n+2}{4} == 3) | (X{n+2}{4} == 5), X{n+1}{4} = 3; else X{n+1}{4} = 6; end
    elseif X{n+1}{2} == X{n+1}{3} - 1
        if (X{n+2}{4} == 1) | (X{n+2}{4} == 3) | (X{n+2}{4} == 5), X{n+1}{4} = 1; else X{n+1}{4} = 2; end
    elseif mod(X{n+1}{2},2) == 1
        if rand < 0.5, X{n+1}{4} = 1; else X{n+1}{4} = 2; end
    else
        if X{n+1}{4} == 1, X{n+1}{4} = 4; else X{n+1}{4} = 5; end
    end
end