function bnode = maximize_params_noisyor(bnode)
% MAXIMIZE_PARAMS_NOISYOR Find the ML parameters of a noisy-or node.
% bnode = maximize_params_noisyor(bnode)

% complete data likelihood = prod_n P(Y=0|X(n))^I(Y(n)=0) * P(Y=1|X(n))^I(Y(n)=1) 
% complete log lik = sum_n (1-Y(n)) log(P(Y=0|X(n)) + Y(n) log(P(Y=1|X(n)) 
% If Y is hidden, taking expectations of this gives
% E[Y(n)] = Pr(Y(n)=1) and E[1-Y(n)] = Pr(Y(n)=0), since Y is a binary indicator variable.
% However, if X is hidden, we can't push E inside the log to get log Pr(Y|E[X])
% Hence it is not correct to replace X with its expected value and use IRLS in this case.
% Instead, we must use EM (see "Learning causal interaction models", Meek and Heckerman, UAI 97).
%
% Alternatively, we can use gradient ascent, but this is numerically difficult.
% The likelihood surface is unimodal in w=log(q), but may have a ridge in q space.
% In both cases, we must use constrained gradient ascent (-inft < w(i) <= 0, or 0 <= q(i) <= 1)
%
% For a Bayesian approach to learning the parameters of a noisy-OR gate, with a Gaussian (!) prior,
% see Diez, "Parameter adjustment in Bayes networks. The generalized noisy-OR gate", UAI 93.

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

mask = logical(bnode.adjustable_mask);

if bnode.observed_parents
  % Pr(Y=0|X=x) = prod_i S(i) ^ X(i) = exp [ sum_i X(i) log S(i) ] = exp [ x' * theta ]
  % where S = inhibit and theta = log(inhibit)
  % The log transformation prevents us from learning a truly deterministic OR (for which inhibit=0).
  theta = log(bnode.inhibit+eps);
  theta = IRLS(bnode.parent_vals, bnode.self_vals, 'noisyor', bnode.max_iter, theta);
  inhibit = exp(theta);
  bnode.inhibit(mask) = inhibit(mask);
else
  bnet = mk_bnet_from_noisyor(bnode.inhibit);
  np = length(bnode.inhibit); % num. parents (including leak)
  N = 2*np+1;
  Cs = 1:np;
  Xs = np+1:2*np;
  E = 2*np+1;

  % We could do all this inside update_ess_noisyor...
  soft_evidence = cell(1, bnode.num_samples);
  for l=1:bnode.num_samples
    ev = cell(1, np+1);
    for i=1:np
      pr = bnode.parent_vals(l,i);
      ev{i} = [1-pr  pr];
    end
    pr = bnode.self_vals(l);
    ev{np+1} = [1-pr  pr];
    soft_evidence{l} = ev;
  end

  if bnode.verbose, fprintf('about to do M step for noisyor node %d ...', bnet.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
    % inhibition prob i = Pr(Ci is on, but Xi is off) 
    if bnode.adjustable_mask(i)
      bnode.inhibit(i) = bnet.bnodes{Xs(i)}.CPT(2,1);
    end
  end
  
end

