function [xsample, tsample, X, P1, randstate] = MEBPsim(X, P1, x, t, n, offspringftn, samplingftn, weightsftn, randstate)

% function [xsample, tsample, X, P1, randstate] = MEBPsim(X, P1, x, t, n, offspringftn, samplingftn, weightsftn, randstate)
%
% simulates an MEBP process
%
% X in is the current state of the process at step 1 (levels 0 to nmax)
% P1 in is the vector of weights for the first crossing at each level 
% x is the position of the process at start of the current (level 0) step
% t is the start time of the current (level 0) crossing
% n is the number of additional steps required
% offspringftn is a function handle for a function that samples from 
%   the offspring distribution p(n)
% samplingftn is a function handle for a function that samples from 
%   the offspring sampling distribution n*p(n)/mu
% 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 n+1
% P1 out is the vector of weights for the first crossing at each level
% xsample is the position of the process at the start of steps 1, ..., n+1
% tsample gives the times the process is at the points in xsample
%
% randstate is an optional parameter for setting the state of the random number generator
%
% offspringftn and samplingftn are functions that take no parameters
%   and return a value in the set {2, 4, 6, ...}
% weightsftn takes a single parameter Z (with value in {2, 4, 6, ...})
%   and returns a vector of length Z of positive reals
%
% row m+1 of X in is (kappa(0,m,1), S(0,m,1), Z(0,m+1,1), alpha(0,m,1), P(0,m+1,1))
% row m+1 of X out is (kappa(0,m,n+1), S(0,m,n+1), Z(0,m+1,n+1), alpha(0,m,n+1), P(0,m+1,n+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 nargin == 9
    rand('state', randstate)
else
    rand('state', sum(100*clock))
end

xsample = zeros(n+1,1);
xsample(1) = x;
tsample = zeros(n+1,1);
tsample(1) = t;

for k = 1:n
	if (X{1}{4} == 1) | (X{1}{4} == 3) | (X{1}{4} == 5) 
        xsample(k+1) = xsample(k) + 1;
    else
        xsample(k+1) = xsample(k) - 1;
    end
    nmax = length(X) - 1;
    for j = 0:nmax, Pt(j+1) = X{j+1}{5}(X{j+1}{2}); end
    tsample(k+1) = tsample(k) + prod(Pt./P1);
	[X, P1] = MEBPexpd(X, P1, samplingftn, weightsftn);
    X = MEBPincr(X, 0, offspringftn, weightsftn);
end

randstate = rand('state');