function CPT = mk_CPT_from_noisyor(q, leak)
% MAKE_CPT_FROM_NOISYOR Compute Pr(X|U1,...,Un) where the Us are the parents (excluding leak)
% CPT = make_CPT_from_noisyor(q, leak)
%
% q(i) is the prob. that the i'th parent will be inhibited (flipped from 1 to 0).
% If leak = 1 (default), q(1) is the leak inhibition probability, and length(q) = n + 1.
% If leak = 0, there is no leak node, and n = length(q).
% The result is CPT(U1,...,Un, X).


if length(q)==1
  CPT = [q  1-q];
  return;
end

if nargin < 2, leak = 1; end

% There are at least 3 different ways of computing this ...
%CPT = f1(q);
CPT = f2(q);
%CPT = f3(q);

n = length(q);
if leak
  CPT = reshape(CPT(2:2:end), twos(1,n)); % skip cases in which the leak is off
else
  CPT = reshape(CPT, twos(1, n));
end

%%%%%%%%

function CPT = f3(q)

n = length(q);
pfail = zeros(n,2);
pfail(:,1) = q(:); % prob. it fails into state 0
CPT = mk_CPT_from_noisy_fn(pfail, 'myor', 2, twos(1,n));

%%%%%

function CPT = f2(q)

n = length(q);
Bn = fliplr(bintab(n)); % all n bit vectors, with the left most column toggling fastest (LSB)
CPT = zeros(2^n, 2);
for i=1:2^n
  u = Bn(i,:);
  CPT(i, 1) = prod(q(find(u)));
  CPT(i, 2) = 1 - CPT(i, 1);
end


%%%%%

function CPT = f1(q)

% Pr(X=0 | U_1 .. U_n) = prod_{i: U_i = on} q_i =  prod_i q_i ^ U_i = exp(u' * log(q_i)
% This method is problematic when q contains zeros
n = length(q);
Bn = fliplr(bintab(n));  % all n bit vectors, with the left most column toggling fastest (LSB)
w = log(q + eps);
M = Bn*w; % M(i) = w . x(i), where x(i) = binary expansion of i
CPT = [exp(M) 1-exp(M)];
