function Gs = simulate_fn(alpha,beta,gamma,maxSize);

zeta = compute_zeta(beta,100000);
%initialize first Graph 
Gs(1).W = 1;
Gs(1).n = 1;
Gs(1).e = 0;

i = 1;

while i <= maxSize
    %Get current graph
    G = Gs(i);
    
    s = sample_starting_node(G,alpha);
    ttl = sample_ttl(G,beta,zeta,G.n-1);
    G_ = [];
    if ttl > G.n-1
        G_ = create_new_node(G,s);
    else
        %run traversal algorithm and return new graph G_
        G_ = traverse(G,s,ttl,gamma);
    end
    i = i+1;
    %add new graph to graph list
    Gs(i) = G_;    
    
end

%**************************************************************************
function G_ = traverse(G,s,ttl,gamma)
visited = zeros(G.n,1);
i = 1;
curr = s;
visited(s) = 1;
path = zeros(ttl+1,1);
path(1) = s;

while i <= ttl
    edges = G.W(curr,:);
    neighbor_idxs = find(edges);
    unvisited_neighbors = find(visited(neighbor_idxs) == 0);
    neighbor_idxs = unvisited_neighbors;
    neighbor_degs = sum(G.W(neighbor_idxs,:),2);
    neighbor_degs = neighbor_degs - 1;
    posdegsIndxs = find(degs > 0);
    if length(posdegsIndxs) == 0        
        break;
    end
    neighbor_degs = neighbor_degs(posdegsIndxs);
    neighbor_idxs = neighbor_idxs(posdegsIndxs);
    sample_probs = neighbor_degs.^gamma;
    sample_probs = sample_probs/sum(sample_probs);
    choose = sample_weighted(sample_probs);
    n_idx = neighbor_idxs(choose);
    visited(n_idx) = 1;
    i = i +1;
    path(i) = n_idx;
    curr = n_idx;
end

if i == ttl+1
    G_ = G;
    G_.e + 1;
    G.W(end,curr) = 1;
    G.W(curr,end) = 1;
else
    G_ = create_new_node(G,curr);
end


%**************************************************************************
function G_ = create_new_node(G,s)
W = spalloc(G.n+1,G.n+1,nnz(G)+1);
W(1:G.n,1:G.n) = G.W;
W(end,s) = 1;
W(s,end) = 1;
G_.W = W;
G_.n = G.n+1;
G_.e = G.e+1;

%**************************************************************************
function s = sample_starting_node(G,alpha)

%out_degs = sum(G.W,1);
%in_degs = sum(G.W,2);
%degs = out_degs + in_degs;
degs = sum(G.W,1);
sample_degs = degs.^alpha;
sample_degs = sample_degs/sum(sample_degs);
sample_degs = cumsum(sample_degs);
s = sample_weighted(sample_degs);

%**************************************************************************
function ttl = sample_ttl(G,beta,zeta,maxSize)
sample = (1:maxSize).^-beta;
sample = [sample, zeta - sum(sample)];
sample = sample/sum(sample);
ttl = sample_weighted(sample);

%**************************************************************************
function choose = sample_weighted(sample)
n = length(sample);
coin = rand;
choose = 1;
for i=1:n
    if sample(i) > coin
        choose = i;
        break;
    end
end