function CPT = mk_CPT_from_noisyfn(pfail, F, child_size)
% MK_CPT_FROM_NOISY_FN Make a noisy version of a deterministic function.
% CPT = mk_CPT_from_noisy_fn(pfail, F, child_size)
% Returns CPT(C1, ..., Cn, E).
%
% We use the algorithm on p211 of Srinivas, "A generalization of the noisy-or model", UAI 93.

if size(pfail,2)==1 % symmetric channel
  pfail = [pfail pfail];
end
np = size(pfail, 1);
parent_size = size(pfail, 2);
psizes = parent_size*ones(1,np);
U = parent_size^np;
CPT = zeros(U, child_size);

XiCPT = cell(1,np);
for i=1:np
  XiCPT{i} = mk_CPT_from_noisy_channel(pfail(i,:), parent_size);
end

if parent_size ~= 2 % general case
  for ui=1:U
    u = ind2subv(psizes, ui);
    for uni=1:U
      un = ind2subv(psizes, uni);
      x = feval(F, un);
      CPT(ui,x) = CPT(ui,x) + prob_noisy2(u, un, XiCPT);
      %CPT(ui, x) = CPT(ui, x) + prob_noisy(u, un, pfail);
    end
  end
else   % Treat binary case separately, for speed
  Bn = fliplr(bintab(np)) + 1; % all n bit vectors, with the left most column toggling fastest (LSB)
  switch F
    case 'myor',
      for ui=1:U
	u = Bn(ui,:);
	for uni=1:U
	  un = Bn(uni,:);
	  x = myor(un);
	  CPT(ui,x) = CPT(ui,x) + prob_noisy2(u, un, XiCPT);
	end
      end
    case 'myand',
      for ui=1:U
	u = Bn(ui,:);
	for uni=1:U
	  un = Bn(uni,:);
	  x = myand(un);
	  CPT(ui,x) = CPT(ui,x) + prob_noisy2(u, un, XiCPT);
	end
      end
    otherwise,
      for ui=1:U
	u = Bn(ui,:);
	for uni=1:U
	  un = Bn(uni,:);
	  x = feval(F, un);
	  CPT(ui,x) = CPT(ui,x) + prob_noisy2(u, un, XiCPT);
	end
      end
  end
end
CPT = reshape(CPT, [psizes child_size]);


%%%%%

function p = prob_noisy2(u, un, XiCPT)
% Returns Pr(X=Un | C=u) for vectors u, un.

np = length(XiCPT);
p = 1;
for i=1:np
  p = p * XiCPT{i}(u(i), un(i));
end

      
%%%%%%%

function p = prob_noisy(u, un, pfail)
% Returns Pr(X=Un | C=u) for vectors u, un.

n = length(u);
P = zeros(1,n);
for i=1:n
  pnofail = 1 - sum(pfail(i,:));
  j = un(i);
  if u(i)==un(i)
    P(i) = pnofail + pfail(i,j);
  else
    P(i) = pfail(i,j);
  end
end
p = prod(P);

