function bnode = maximize_params_noisyfn(bnode)
% MAXIMIZE_PARAMS_NOISYFN Find the ML parameters of a noisy-fn node.
% bnode = maximize_params_noisyfn(bnode)
%
% We follow the technique discussed in
% "Learning causal interaction models", Meek and Heckerman, UAI 97.
% See the README file for details.

% The CPT of node Xi looks like this, where p(j) = pfail(i,j),
% and . means 1 - the sum of the other cols (child vals)
% C\X  1     2   3
%  1   .   p(2) p(3)
%  2  p(1) .    p(3)
%  3  p(1) p(2) .
% Hence we can estimate pfail(i,j) by finding a p(j) in this table.
% For non-binary nodes, this is not optimal, since there is more than one such p(j):
% we should have pooled the stats of each copy, since they are tied.
% If this is statistical inefficiency is a problem, you should create a noisy_node type to model
% a noisy channel, and use it instead of tabular nodes for the Xi's.
% (A noisy_node behaves just like a noisyfn_node where F is the identity fn.)

% We treat symmetric channels similarly, i.e., non-optimally.
% e.g., for a binary symmetric channel, the CPT for Xi looks like this
% 
% C\X  0   1
%  0   1-p p
%  1   p   1-p
% 
% where p = pfail(i,1) = pfail(i,2).
% We estimate p by examining cell (0,1), even though there are two copies.
% (Tables respect the constraint that the rows sum to 1.)

% We treat clamped (non-adjustable) parameters in a rather cheesy way:
% Inside the EM routine, they are unconstrained, but afterwards, we set them to
% their original values.


bnet = mk_bnet_from_noisyfn(bnode.pfail, bnode.fn, bnode.self_size, bnode.prior);
np = length(bnode.parents); % num. parents
Cs = 1:np;
Xs = np+1:2*np;
E = 2*np+1;

% We could do this in update_ess_noisyfn...
N = bnode.num_samples;
soft_evidence = cell(1, N);
for l=1:N
  % each row gets its own cell
  ev = num2cell(bnode.parent_distrib(:,:,l), 2);
  ev{end+1} = bnode.self_distrib(:,l)'; % might be a different length to the parent distrib
  soft_evidence{l} = ev;
end

if bnode.verbose, fprintf('about to do M step for noisyfn node %d\n', bnode.self); end
bnet = learn_params(bnet, [], [], bnode.max_iter, bnode.thresh, bnode.verbose, soft_evidence, [Cs E]);
if bnode.verbose, fprintf('done\n'); end

for i=1:np
  if bnode.symmetric
    if bnode.adjustable_mask(i)
      bnode.pfail(i) = bnet.bnodes{Xs(i)}.CPT(1,2);
    end
  else
    pf = zeros(1, bnode.parent_size);
    pf(1) = bnet.bnodes{Xs(i)}.CPT(2,1);
    pf(2:end) = bnet.bnodes{Xs(i)}.CPT(1,2:end);
    mask = logical(bnode.adjustable_mask(i,:));
    bnode.pfail(i,mask) = pf(mask);
  end
end

