function [bnet, ll] = enter_evidence(bnet, vals, obs_nodes, soft_vals, soft_nodes)
% ENTER_EVIDENCE Clear old evidence, add new evidence, and propagate its effects.
%
% [BNET, LL] = ENTER_EVIDENCE(BNET, VALS) adds the the evidence in VALS and calls an inference
% engine. If BNET is a static BN, VALS contains the values of the observed nodes concatenated in order.
% If BNET is a DBN, column t of VALS contains the values for timeslice t.
% Note that, for discrete nodes, the values must be positive integers (0s are disallowed).
%
% [BNET, LL] = ENTER_EVIDENCE(BNET, VALS, OBS_NODES) specifies what nodes are observed.
% (OBS_NODES defaults to the nodes that were declared to be observed when BNET was made.)
% If BNET is a DBN, and OBS_NODES is a vector, it is assumed that the same nodes are observed in every timeslice.
% If different nodes are observed in every time slice (e.g., because of missing information)
% VALS and OBS_NODES should be cell arrays: OBS_NODES{t} specfies which nodes are
% observed in slice t, and VALS{t} is a vector of their values.
% (This is only possible if bnet.algo == 'bk' or 'jtree'.)
%
% [BNET, LL] = ENTER_EVIDENCE(BNET, VALS, OBS_NODES, SOFT_VALS, SOFT_NODES) provides soft
% (virtual) evidence on the specified nodes. (A node is said to have soft  evidence if
% we don't know its exact value, but have a prob. distrib. over its values.)
% SOFT_NODES specfies which nodes have soft evidence, and SOFT_VALS{I} is a prob. distribution
% over the values for node SOFT_NODES(I). Currently only discrete nodes can have soft evidence,
% and you must use the jtree algorithm. (Soft evidence is used in maximize_params_noisyor.)
%
% After the evidence has been propagated, marginals on nodes can be computed
% using marginal_nodes/ marginal_family.
%
% Each inference algorithm takes various optional arguments.
% These are passed in as fields, as explained below. They should be set before calling this function.
% 
% - jtree
%     bnet.maximize = 1 if we want to find the most probable assignment (c.f. Viterbi path in an HMM).
%     i.e., replace the sum operator with a max operator (Default: maximize = 0).
%     bnet.mpa(i) contains the most probable value of node i (in the unrolled network).
%     Remember to reset bnet.maximize = 0 after calling.
%
% - pearl
%     When the graph has loops (undirected cycles), message passing will continue until the
%     change in beliefs falls below pearl.threshold, or pearl.max_iter iterations are exceeded,
%     whichever occurs first. The number of iterations is returned in pearl.niterations.
%     bnet.name contains the name of a file to which the beliefs are written after each iteration.
%     bnet.maximize = 1 is as for jtree.
%
% - cutset
%     cutset.cutset_nodes should contain a valid loop cutset
%
% - lw, ls
%     sample.num_samples determines the number of samples to use.
%
% - sof
%     sof.num_samples determines the number of samples to use.
%     sof.active = 1 means do active sampling (default: always use exact num_samples samples)
%
% - bk, hmm, lds
%     bnet.filter_only = 1 means don't run the backwards step. Default: do full smoothing.


raw_vals = vals; % before we convert to cell array form (used by HMM/LDS algos)

if ~iscell(vals)
  if bnet.is_dbn
    vals = num2cell(vals, 1); % each column gets its own cell
  else
    vals = {vals(:)};
  end
end

if bnet.is_dbn
  T = length(vals);
  if exist('obs_nodes') & iscell(obs_nodes)
    T = max(T, length(obs_nodes));
  end
  if bnet.num_slices ~= T
    bnet = unroll_dbn(bnet, T);
    if strcmp(bnet.algo, 'jtree')
      % since the size (length) of the network has changed
      bnet.jtree = mk_jtree(bnet);
    end
  end
else
  T = 1;
end

if ~exist('obs_nodes')
  if bnet.is_dbn
    obs_nodes = bnet.observed_slice;
  else
    obs_nodes = bnet.observed_nodes(:);
  end
end

if ~iscell(obs_nodes)
  if bnet.is_dbn
    % each column is the same, and each column gets its own cell
    obs_nodes = num2cell(repmat(obs_nodes(:), [1 T]), 1); 
  else
    obs_nodes = {obs_nodes(:)};
  end
end

if ~exist('soft_nodes'),  soft_nodes = []; end
if ~exist('soft_vals'),  soft_vals = {}; end


num_nodes = length(bnet.adj_mat);
all_nodes = 1:num_nodes;
ss = bnet.slice_size;

% Change to an unrolled numbering system
old_obs_nodes = bnet.observed_nodes;
bnet.observed_nodes = [];
for t=1:T
  bnet.observed_nodes = [bnet.observed_nodes obs_nodes{t}(:)' + (t-1)*ss];
end

% ~isequal(empty 1x0 matrix, [])
if prod(size(bnet.observed_nodes))==0, bnet.observed_nodes = []; end
bnet.hidden_nodes = setdiff(all_nodes, bnet.observed_nodes);
if prod(size(bnet.hidden_nodes))==0, bnet.hidden_nodes = []; end

% It is not sufficient just to pass in vals to the inference algorithm,
% since add_evidence_to_marginal also needs it.
% So we must store it inside bnet as a concatentation of all the columns.
if isempty(vals)
  all_vals = [];
else
  %all_vals = cat(1, vals{:}); % fails if some slices have differing amts of data
  all_vals = [];
  for i=1:length(vals)
    all_vals = [all_vals vals{i}];
  end
end
bnet.values = blockmatrix(1, [], 1, bnet.blocks(:)); % clear old evidence
%%%bnet.values = blockmatrix(1, bnet.blocks(:)); % clear old evidence
bnet.values(bnet.observed_nodes) = all_vals(:)'; % add new evidence

nnodes = length(bnet.adj_mat);
bnet.obs_mask = zeros(1, nnodes);
bnet.obs_mask(bnet.observed_nodes) = 1; % bit vectors are faster than sets

% Constraint propagation.
% Currently we implement a very limited form, where deterministic children become observed iff
% their immediate parents are observed.
for i=1:nnodes
  [determined, val] = determine_value_node(bnet, i);
  if determined
    bnet.values(i) = val;
    bnet.obs_mask(i) = 1;
    % the size of observed nodes will be temporarily changed inside mk_potential
  end
end
bnet.observed_nodes = find(bnet.obs_mask);

% If we wait until the evidence has arrived before making the jtree,
% we can disconnect parents to multiplexer nodes if their gate values are known.
% This can result in drastically smaller cliques.
if strcmp(bnet.algo, 'jtree')
  if bnet.wait_till_evidence
    assert(~bnet.is_dbn);
    % Make a copy of the current graph before we modify its structure
    original_dag = bnet.adj_mat;
    for i=1:nnodes
      [changed, new_ps] = change_parents_node(bnet, i);
      if changed
	old_ps = parents(bnet.adj_mat, i);
	bnet.adj_mat(old_ps, i) = 0;
	bnet.adj_mat(new_ps, i) = 1;
	% this may invalidate any bnode.parents fields...
	% (Fortunately, we don't store bnode.children!)
      end		
    end
    bnet.new_adj_mat = bnet.adj_mat;
    bnet.jtree = mk_jtree(bnet);
  end
end

%if strcmp(bnet.algo, 'jtree') & bnet.is_dbn & (bnet.num_slices ~= T)
%  % since the size (length) of the network has changed
%  bnet.jtree = mk_jtree(bnet);
%end



% Is the hidden state all discrete or all Gaussian or mixed?
if strcmp(bnet.algo, 'jtree') |  strcmp(bnet.algo, 'bk')
  if isequal(bnet.observed_nodes(:), all_nodes(:))
    bnet.jtree.pot_type = 'd'; % could also use 'g' or 'cg'
  else
    if subset(bnet.hidden_nodes, bnet.discrete_nodes) % [] is a subset of everything
      bnet.jtree.pot_type = 'd';
    elseif subset(bnet.hidden_nodes, bnet.cts_nodes);
      bnet.jtree.pot_type = 'g';
    else
      %disp('network has discrete and cts hidden state');
      bnet.jtree.pot_type = 'cg';
    end
  end
  if bnet.maximize
    assert(bnet.jtree.pot_type == 'd');
    bnet.jtree.pot_type = 'm'; % max tables
  end
end
if strcmp(bnet.algo, 'bk')
  bnet.jtree1.pot_type = bnet.jtree.pot_type;
end


% Are there any links from hidden cts to discrete (logistic) nodes?
% If so, we must use a variational algo (even if the discrete child is observed).
% See "A Variational Approximation for Bayesian networks with discrete and continuous
% latent variables", K. Murphy, UAI99.

% For the latest details on hybrid inference (discrete parents, cts children), see
% S. L. Lauritzen and F. Jensen. Stable local computation with conditional Gaussian distributions. 
% Research Report R-99-2014. Department of Mathematical Sciences, Aalborg University, 1999
% http://www.math.auc.dk/~steffen/publications.html

if strcmp(bnet.algo, 'jtree')   
  var = 0;
  if bnet.jtree.pot_type ~= 'd'
    for i=1:num_nodes
      if ismember(i, bnet.cts_nodes) & ismember(i, bnet.hidden_nodes)
	cs = children(bnet.adj_mat, i);
	if ~isempty(intersect(cs, bnet.discrete_nodes)) % overlap(cs, bnet.discrete_nodes)
	  var = 1;
	  break;
	end
      end
    end
  end
end


switch bnet.algo
  case 'jtree'
    if var
      warning('the code to handle discrete children with hidden cts parents is very experimental!');
      [bnet, ll] = propagate_var_jtree(bnet);
    else
      [bnet, ll] = propagate_jtree(bnet, soft_nodes, soft_vals);
    end
  case 'hmm'
    [bnet, ll] = propagate_hmm(bnet, raw_vals);
  case 'lds'
    [bnet, ll] = propagate_lds(bnet, raw_vals);
  case 'bk'
    % propagate_dbn does not read or write bnet.values or bnet.observed_nodes
    [bnet, ll] = propagate_dbn(bnet, obs_nodes, vals);
  case 'pearl'
    [bnet, ll] = propagate_pearl(bnet);
  case 'cutset'
    [bnet, ll] = propagate_cutset(bnet);
  case 'lw'
    [bnet, ll] = likelihood_weighting(bnet);
  case 'ls',
    [bnet, ll] = logic_sampling(bnet);
  case 'sof'
    [bnet, ll] = sof(bnet, obs_nodes, vals);
  otherwise,
    error(['unrecognized inference algorithm ' bnet.algo]);
end

if strcmp(bnet.algo, 'jtree') & bnet.wait_till_evidence
  % restore to original structure so that marginal_family and future use will work
  bnet.adj_mat = original_dag;
end
