function [FB,EB] = pbnBayesPreds(v,k,F,bi)

% [FB,EB] = pbnBayesPreds(v,k,F) - Bayes steady-state predictors for a PBN
%
% Function computes the Bayes steady-state prediction error for all the 
% variables and all the predictor variable combinations. Note that a 
% variable is not allowed to be used as a predictor variable for itself. In
% case of tie, a randomly selected function with the minimum Bayes error is
% returned.
%
% INPUT:
% v     - The stationary distribution of a PBN, i.e., f(1) = 
%         Pr{X = 00...00}, f(2) = Pr{X = 00...01}, ..., f(2^n) = 
%         Pr{X = 11...11}.
% k     - The number of variables used in each predictor function.
% F     - The set of Boolean predictors to be used in the inference. If F 
%         is the empty matrix, then the function class is considered to 
%         contain all k-variable Boolean functions. F is either a
%         (2^k)-by-nf binary matrix or 1-by-nf row vector of integers (see
%         also the description of the next argument bs), where k is the 
%         number of variables in each function and nf is the number of 
%         functions.
%         1.    The case of (2^k)-by-nf binary matrix: Let f = F(:,j) be
%               the j:th column of F (i.e., the j:th truth table in F. 
%               Then, f(0) defines the output value for the input vector
%               00...00, f(1) for the input 00...01, f(2) for the input 
%               00...10, ..., and f(2^k) for the input 11...11. Input 
%               vectors are interpreted such that the left most bit defines
%               the value of the first input variable, the second bit from 
%               the left defines the value of the second input variable, 
%               ..., and the right most bit defines the value of the last 
%               (k:th) input variable.
%         2.    The case of 1-by-nf row vector of integers: Let f = F(j) be
%               the j:th element of F. Then, the first bit (as obtained by 
%               the bitget command bitget(f,1)) defines the output value 
%               for the input vector 00...00, the second bit bitget(f,2) 
%               defines the output for the input 00...01, ..., and the 
%               2^k:th bit bitget(f,2^k) defines the output for the input
%               11...11. Thus, the regular truth table presentation of the
%               j:th function can be obtained by f = bitget(F(j),[1:2^k])'.
%               Note that only the cases k<=5 can be handled by this
%               convention.
% bi    - A bit (0/1) indicating that whether the Boolean functions in F 
%         are represented in the form of standard binary truth tables (0) 
%         or (encoded) integers (1). (This can be used to distinguish between 
%         constant functions and the integer presentations).
% OUTPUT:
% FB    - A 3-D matrix of the Bayes predictors for the all variables and 
%         all the predictor variable combinations. FB has size
%         (2^k)-by-nchoosek(n-1,k)-by-n, where n is the number of variables
%         in the PBN. FB(:,:,i) defines the best functions for the i:th 
%         node. Each column in FB, i.e., FB(:,i,j) is interpreted as the 
%         columns in F (see above).
% EB    - A 2-D matrix of the Bayes errors for all the variables and all 
%         the predictor variable combinations. EB has size 
%         nchoosek(n-1,k)-by-n.

% Functions used: margpdf.m, randintex.m

% 26.08.2005 by Harri L鋒desm鋕i. Modified from pbnBayesPred.m


%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Define and initialize some variables.
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

n = round(log2(length(v))); % The number of nodes.
kk = 2^k; % Needed often.

combnum = nchoosek(n-1,k); % The number of different variable combinations.

% All variable combinations will be generated in advance. Check that he
% number of combinations is "samll" enough. This will work only for small
% enough PBNs.
if combnum>20000 % Limit the number of possible combinations.
    error('Too many variable combinations...')
end % if combnum>20000

starti = 1;
stopi = combnum;

% Initialize the output matrices.
FB = zeros(kk,combnum,n);
EB = zeros(combnum,n);


%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% The main loop separately for unconstrained and constrained case.
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% If F is an empty matrix, then the function class is considered to contain
% all k-variables Boolean functions (unconstrained case).
if isempty(F)
    
    % Run through all the nodes.
    for i=1:n
        
        % Generate all variable combinations in advance. Remove the current
        % node (target node) from the set of predictors.
        IAll = nchoosek([1:i-1,i+1:n],k);
        
        %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        % Run through all variable combinations.
        %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        for j=starti:stopi
            
            % The current variable combinations (in lexicographical ordering).
            I = IAll(j,:);
            
            % Compute the corresponding marginal distribution (including
            % the target node itself).
            vv = margpdf(v,[i,I]);
            
            % Find the optimal (Bayes) k-variable predictor for the current
            % input variable combination.
            [OptErr,OptF] = min([vv(1:kk);vv(kk+1:2*kk)]);
            OptF = 2 - OptF;
            % All output bits having tie are set uniformly randomly.
            Ties = (vv(:,1:kk)==vv(:,kk+1:2*kk));
            OptF(Ties) = (rand(1,sum(Ties(:))))>0.5;
            % The corresponding probability of error.
            OptErr = sum(OptErr);
            
            FB(:,j,i) = OptF';
            EB(j,i) = OptErr;
            
        end % for i=starti:stopi
    end % for i=1:n
    
else
    
    % If F is non-empty matrix, then the function class is defined by F.
    
    nf = size(F,2); % The number of functions in F.
    Oneskknf = ones(kk,nf);
    Zeroskknf = zeros(kk,nf);
    Ones1nf = ones(1,nf);
    
    % Run through all the nodes.
    for i=1:n
        
        minerr = realmax; % Keep track of the minimum error.
        
        % Generate all variable combinations in advance. Remove the current
        % node from the set of predictors.
        IAll = nchoosek([1:i-1,i+1:n],k);
        
        %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        % Run through all variable combinations.
        %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        for j=starti:stopi
            
            % The current variable combinations (in lexicographical ordering).
            I = IAll(j,:);
            
            % Compute the corresponding marginal distribution (including
            % the target node itself).
            vv = margpdf(v,[i,I]);
            
            % Compute the error of all the functions.
            Err = sum(F.*(vv(1:kk)'*Ones1nf)) + ...
                sum((1-F).*(vv(kk+1:2*kk)'*Ones1nf));
            
            % Find the optimal (Bayes) k-variable predictor for the current
            % input variable combination.
            OptErr = min(Err);
            
            if OptErr<minerr
                % Select one of the optimal functions randomly.
                ind = find(Err==OptErr);
                ind = ind(randintex(1,1,length(ind)));
                FB(:,i) = F(:,ind);
                varF(:,i) = I';
                EB(i) = OptErr;
                minerr = OptErr;
            end % if OptErr<minerr
                        
        end % for i=starti:stopi
    end % for i=1:n
end % if isempty(F)