function pst = build_pst(data, model_score)

nfval = 3;

pst{1}.x = [];
incomplete = [1];
leafs = zeros([nfval 1]);

score = score_pst(pst, data, model_score);

leaf_added = 1;

while(leaf_added)

    K = length(pst);

    leaf_added = 0;
    best_score = score;
    best_node = 0;
    best_edge = 0;
    
    for i = 1:length(incomplete)
        edges = find(leafs(:,i) == 0);

        % force adding the last edge of each incomplete node
        if(length(edges) == 1)
            leaf_added = 1;
            best_node = i;
            best_edge = edges(1);
            best_score = score;
            break;
        end

        for j = 1:length(edges)
            new_pst = pst;
            new_pst{K+1}.x = [pst{incomplete(i)}.x edges(j)-1];
            new_score = score_pst(new_pst, data, model_score);
            if(new_score > best_score)
                best_score = new_score;
                best_node = i;
                best_edge = edges(j);
                leaf_added = 1;
            end
        end
    end
    
    % add the best leaf
    if leaf_added
        pst{K+1}.x = [pst{incomplete(best_node)}.x best_edge-1];
        score = best_score;
        leafs(best_edge, best_node) = 1;
        incomplete = [incomplete K+1];
        leafs = [leafs zeros([nfval 1])];

        % remove the best node if completed
        if(sum(leafs(:,best_node)) == nfval)
            [incomplete, I] = setdiff(incomplete, incomplete(best_node));
            leafs = leafs(:,I);
        end
    end
end
