function pot = add_evidence_to_marginal(bnet, m)
% ADD_EVIDENCE_TO_MARGINAL Add the values of observed nodes to a potential.
% pot = add_evidence_to_marginal(bnet, m)
%
% We proceed as in the following example.
% Suppose m has domain [1,2,3,4] where 1 and 3 are discrete, 2 and 4 are cts,
%
% Case 1: m is discrete (because all the hidden nodes are discrete)
% 1 is discrete and is observed to have value 5.
% 2 is cts and is observed to have value [5 42]
% 3 is discrete and is hidden (but has 3 possible values)
% 4 is cts and is observed to have value [1].
%
% Hence m.sizes = [1 1 3 1]
% index0 = {:,1,:,1} to extract the discrete nodes 1 and 3.
% index1 = {5,:} to insert a position compatible with the discrete evidence
% pot.T(5,:) = m.T(1,1,:,1)
% pot.mu(:,5,:) = observed values of [X2 X4] = cobs_vals
% pot.Sigma = zeros since all cts nodes are observed
%
% Case 2: m is cts (because all the hidden nodes are cts)
% 1 is discrete and is observed to have value 5.
% 2 is cts and is observed to have value [5 42]
% 3 is discrete and is observed to have value 1.
% 4 is cts and is hidden (but is a vector of length 3)
%
% Hence m.sizes = [0 0 0 3]
% indexo = [1 2] = indices containing X2 (cts_obs_dom)
% indexh = [3 4 5] = indices containing X4 (cts_hid_dom)


if ~any(bnet.obs_mask(m.domain)) % isempty(odom)
  pot = m;
  if m.type=='g' | strcmp(m.type, 'cg')
    pot.mu = m.mu;
    pot.Sigma = m.Sigma;
  end
  return;
end

odom = intersect(m.domain, bnet.observed_nodes);
ddom = intersect(m.domain, bnet.discrete_nodes); 
cdom = intersect(m.domain, bnet.cts_nodes);
hdom = intersect(m.domain, bnet.hidden_nodes);
dodom = intersect(odom, bnet.discrete_nodes);
codom = intersect(odom, bnet.cts_nodes);
dhdom = intersect(hdom, bnet.discrete_nodes);
chdom = intersect(hdom, bnet.cts_nodes);
dvals = bnet.values(dodom);
cvals = bnet.values(codom);
index1 = mk_multi_index(length(ddom), find_equiv_posns(dodom, ddom), dvals);
ns = bnet.node_sizes;
dsizes = ns(ddom);
dsize = prod(ns(ddom));
csize = sum(ns(cdom));

pot.domain = m.domain;
pot.type = 'cg';
pot.sizes = ns(m.domain);
if ~isempty(ddom)
  pot.T = myzeros(dsizes);
else
  pot.type = 'g';
end
if ~isempty(cdom)
  pot.mu = myzeros([csize dsizes]);
  pot.Sigma = myzeros([csize csize dsizes]);
else
  pot.type = 'd';
end


switch m.type
  case {'d', 'm'}, % case 1
    index0 = mk_multi_index(length(m.domain), find_equiv_posns(cdom, m.domain), ones(1, length(cdom)));
    pot.T(index1{:}) = m.T(index0{:});
    ns(bnet.observed_nodes) = 1;
    if ~isempty(cdom)
      pot.mu(:, index1{:}) = repmat(full(cvals(:)), [1 ns(ddom)]);
    end
  case 'g', % case 2
    assert(m.subtype == 'm');
    [mu, Sigma] = add_ev_to_gaussian(m.mu, m.Sigma, chdom, codom, cdom, cvals, ns(cdom));
    % Insert into discrete table (if there is one)
    if isempty(ddom)
      pot.mu = mu;
      pot.Sigma = Sigma;
    else
      pot.T(index1{:}) = 1; % Pr(X=x|E)=1 if X=x is part of the evidence
      pot.mu(:,index1{:}) = mu;
      pot.Sigma(:,:,index1{:}) = Sigma;
    end
  case 'cg', % case 3
    mask = find_equiv_posns(dodom, ddom);
    for i=1:dsize
      sub = ind2subv(dsizes, i);
      sub(mask) = 1;
      j = subv2ind(m.dsizes, sub);
      if length(cvals)>0
	[mu, Sigma] = add_ev_to_gaussian(m.pot{j}.mu, m.pot{j}.Sigma, chdom, codom, cdom, cvals, ns(cdom));
      else
	mu = m.pot{j}.mu;
	Sigma = m.pot{j}.Sigma;
      end
      % Insert into discrete table (if there is one)
      if isempty(ddom)
	pot.mu = mu;
	pot.Sigma = Sigma;
      else
	pot.T(i) = exp(m.pot{j}.logp);
	%pot.T(i) = pot.T(i) / normal_coef(m.pot{j}.Sigma)  % COEF
	pot.mu(:,i) = mu;
	pot.Sigma(:,:,i) = Sigma;
      end
    end
  otherwise,
    error(['bad type ' m.type]);
end


  
