wolffd@0: function [engine, loglik] = enter_evidence(engine, evidence, varargin)
wolffd@0: % ENTER_EVIDENCE Add the specified evidence to the network (cond_gauss)
wolffd@0: % [engine, loglik] = enter_evidence(engine, evidence, ...)
wolffd@0: %
wolffd@0: % evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: ns = bnet.node_sizes(:);
wolffd@0: observed = ~isemptycell(evidence);
wolffd@0: onodes = find(observed);
wolffd@0: hnodes = find(isemptycell(evidence));
wolffd@0: engine.evidence = evidence;
wolffd@0: 
wolffd@0: % check there are no C->D links where C is hidden
wolffd@0: pot_type = determine_pot_type(bnet, onodes);
wolffd@0: 
wolffd@0: dhid = myintersect(hnodes, bnet.dnodes);
wolffd@0: S = prod(ns(dhid));
wolffd@0: T = zeros(S,1);
wolffd@0: 
wolffd@0: N = length(bnet.dag);
wolffd@0: mu = cell(1,N);
wolffd@0: Sigma = cell(1,N); 
wolffd@0: cobs = myintersect(bnet.cnodes, onodes);
wolffd@0: chid = myintersect(bnet.cnodes, hnodes);
wolffd@0: ens = ns;
wolffd@0: ens(cobs) = 0;
wolffd@0: for j=chid(:)'
wolffd@0:   mu{j} = zeros(ens(j), S);
wolffd@0:   Sigma{j} = zeros(ens(j), ens(j), S);
wolffd@0: end
wolffd@0:  
wolffd@0: for i=1:S
wolffd@0:   dvals = ind2subv(ns(dhid), i);
wolffd@0:   evidence(dhid) = num2cell(dvals);
wolffd@0:   [sub_engine, loglik] = enter_evidence(engine.sub_engine, evidence);
wolffd@0:   for j=chid(:)'
wolffd@0:     m = marginal_nodes(sub_engine, j);
wolffd@0:     mu{j}(:,i) = m.mu;
wolffd@0:     Sigma{j}(:,:,i) = m.Sigma;
wolffd@0:   end
wolffd@0:   T(i) = exp(loglik);
wolffd@0: end
wolffd@0: 
wolffd@0: [T, lik] = normalise(T);
wolffd@0: loglik = log(lik);
wolffd@0: 
wolffd@0: engine.T = T;
wolffd@0: engine.mu = mu;
wolffd@0: engine.Sigma = Sigma;
wolffd@0: 
wolffd@0: dnodes = bnet.dnodes;
wolffd@0: dobs = myintersect(dnodes, onodes);
wolffd@0: ens(dobs) = 1;
wolffd@0: engine.joint_dmarginal = dpot(dnodes, ens(dnodes), myreshape(engine.T, ens(dnodes)));
wolffd@0: 
wolffd@0: engine.onodes = onodes;