function [xsample, tsample] = MEBP(offspring_dist, offspring_param, weights_dist, weights_param, crossing_dist, crossing_param, n)

% front end for function MEBPsim
%
% requires anonymous functions from Matlab version 7 (R14)
%
% n is number of sample points required
% xsample and tsample give position and time of simulated process
% (both starting at 0 and length n + 1)
%
% offspring_dist is a text string indicating which distribution to use for
% the offspring; offspring_param are the parameters of that distribution
% current options are 
%   'geo' for 2*geometric(p) + 2
%     param: p
%   'geo2' for 2*X where X is the time to success for a sequence of trials
%     where the first success probability is a and subsequently is b
%     param: [a, b]
%   'twovalue' for a r.v. that takes on values 2 and x and has mean 2^(1/H)
%     param: [x, H]
%
% weights_dist is a text string indicating which distribution to use for
% the weights; weights_param are the parameters of that distribution
% current options are
%   'const' for constant weights
%     param: []
%   'indgam' for independent gamma(a, 1) weights
%     param: a
%   'indlogn' for independent lognormal(mu, si2) weights
%     param: [mu, si2]
%
% crossing_dist is a text string indicating which distribution to use for
% the crossings; crossing_param are the parameters of that distribution
% current option is
%   'indpt' for crossings independent of the size of the parent crossing
%   param: p, up-down with proba p and down-up with proba 1-p
%    
%
% examples
%   [x, t] = MEBP('geo', .5, 'const', [], 'indpt', .5, 10000); % simple random walk
%   [x, t] = MEBP('geo', .7, 'const', [], 'indpt', .5, 10000); % long range dept
%   [x, t] = MEBP('geo', .5, 'indgam', 2, 'indpt', .5, 10000); % BM with multifractal
%     time change
%   [x, t] = MEBP('geo', .5, 'indlogn', [0 1], 'indpt', .5, 10000); % BM with
%     multifractal time change
%
% Owen Jones 17 January 2007, 18 June 2007, Geoffrey Decrouez 25 July 2008

if strcmp(offspring_dist, 'geo')
    offspringftn = @(x)2*geornd(offspring_param)+2;
    samplingftn = @(x)geosample(offspring_param);
elseif strcmp(offspring_dist, 'geo2')
    offspringftn = @(x)geo2rnd(offspring_param(1), offspring_param(2));
    samplingftn = @(x)geo2sample(offspring_param(1), offspring_param(2));
elseif strcmp(offspring_dist, 'twovalue')
    offspringftn = @(x)twovaluernd(offspring_param(1), offspring_param(2));
    samplingftn = @(x)twovaluesample(offspring_param(1), offspring_param(2));
end

if strcmp(weights_dist, 'const')
    weightsftn = @(Z)ones(1,Z);
elseif strcmp(weights_dist, 'indgam')
    weightsftn = @(Z)gamrnd(weights_param, 1, [1, Z]);
elseif strcmp(weights_dist, 'indlogn')
    weightsftn = @(Z)lognrnd(weights_param(1), weights_param(2), [1, Z]); 
end

if strcmp(crossing_dist, 'indpt')
    alphaftn = @(Z)crossrnd(crossing_param,Z);
    alphaftn2 = @(Z,Ztype)crossrnd2(crossing_param,Z,Ztype);
    alphacondftn = @(Z,S,Stype)crosscondrnd(crossing_param,Z,S,Stype);
end

[X, P1] = MEBPinit(samplingftn, weightsftn, alphaftn);
[xsample, tsample, X, P1] = MEBPsim(X, P1, 0, 0, n, offspringftn, samplingftn, weightsftn, alphaftn2, alphacondftn);


function n = geosample(p)
r = rand;
n = 1;
prb = p^2;
while r > prb
    n = n + 1;
    prb = prb + n*(1-p)^(n-1)*p^2;
end
n = 2*n;

function n = geo2rnd(a, b)
r = rand;
n = 1;
prb = a;
while r > prb
    n = n + 1;
    prb = prb + (1-b)^(n-2)*(1-a)*b;
end
n = 2*n;

function n = geo2sample(a, b)
r = rand;
n = 1;
m = (1 - a + b)/b;
prb = a/m;
while r > prb
    n = n + 1;
    prb = prb + n*(1-b)^(n-2)*(1-a)*b/m;
end
n = 2*n;

function n = twovaluernd(x, H)
mu = 2^(1/H);
p = (mu-2)/(x-2);
if rand < 1-p, n = 2; else n = x; end

function n = twovaluesample(x, H)
mu = 2^(1/H);
p = (mu-2)/(x-2);
if rand < 2*(1-p)/mu, n = 2; else n = x; end

function alpha=crossrnd(crossing_param,Z)
if rand<crossing_param
    Ztype= 'u';
else
    Ztype='d';
end
alpha= zeros(1,Z);
if Z>2
    for i=1:((Z/2)-1)
        if rand <crossing_param
            alpha(2*i-1)=1;
            alpha(2*i)=4;
        else
            alpha(2*i-1)=2;
            alpha(2*i)=5;
        end
    end
end
if Ztype=='u'
    alpha(Z-1)=1;
    alpha(Z)=3;
else
    alpha(Z-1)=2;
    alpha(Z)=6;
end

function alpha=crossrnd2(crossing_param,Z,Ztype)
if Z>2
    for i=1:(Z/2-1)
        if rand <crossing_param
            alpha(2*i-1)=1;
            alpha(2*i)=4;
        else
            alpha(2*i-1)=2;
            alpha(2*i)=5;
        end
    end
end
if Ztype=='u'
    alpha(Z-1)=1;
    alpha(Z)=3;
else
    alpha(Z-1)=2;
    alpha(Z)=6;
end

function alpha=crosscondrnd(crossing_param,Z,S,Stype)
alpha= crossrnd(crossing_param,Z);
if Stype=='u'
    if (S==Z-1)|(S==Z)
        alpha(Z-1)=1;
        alpha(Z)=3;
    elseif mod(S,2)==1
        alpha(S)=1;
        alpha(S)=4;
    else
        alpha(S-1)=2;
        alpha(S)=5;
    end
else
    if (S==Z-1)|(S==Z)
        alpha(Z-1)=2;
        alpha(Z)=6;
    elseif mod(S,2)==1
        alpha(S)=2;
        alpha(S)=5;
    else
        alpha(S-1)=1;
        alpha(S)=4;
    end    
end