wolffd@0: function [engine, ll, niter] = enter_evidence(engine, evidence, varargin)
wolffd@0: % ENTER_EVIDENCE Propagate evidence using belief propagation
wolffd@0: % [engine, ll, niter] = enter_evidence(engine, evidence, ...)
wolffd@0: %
wolffd@0: % The log-likelihood is not computed; ll = 0.
wolffd@0: % niter contains the number of iterations used (if engine.protocol = 'parallel')
wolffd@0: %
wolffd@0: % The following optional arguments can be specified in the form of name/value pairs:
wolffd@0: % [default value in brackets]
wolffd@0: %
wolffd@0: % maximize - 1 means use max-product, 0 means use sum-product [0]
wolffd@0: % exclude  - list of nodes whose potential will not be included in the joint [ [] ]
wolffd@0: %
wolffd@0: % e.g., engine = enter_evidence(engine, ev, 'maximize', 1)
wolffd@0: 
wolffd@0: ll = 0;
wolffd@0: exclude = [];
wolffd@0: maximize = 0;
wolffd@0: 
wolffd@0: if nargin >= 3
wolffd@0:   args = varargin;
wolffd@0:   nargs = length(args);
wolffd@0:   for i=1:2:nargs
wolffd@0:     switch args{i},
wolffd@0:      case 'exclude', exclude = args{i+1};
wolffd@0:      case 'maximize', maximize = args{i+1};
wolffd@0:      otherwise,
wolffd@0:       error(['invalid argument name ' args{i}]);
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: engine.maximize = maximize;
wolffd@0: 
wolffd@0: if ~isempty(engine.filename)
wolffd@0:   engine.fid = fopen(engine.filename, 'w');
wolffd@0:   if engine.fid == 0
wolffd@0:     error(['can''t open ' engine.filename]);
wolffd@0:   end
wolffd@0: else
wolffd@0:   engine.fid = [];
wolffd@0: end
wolffd@0: 
wolffd@0: gdl = engine.gdl;
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: 
wolffd@0: ndoms = length(gdl.doms);
wolffd@0: ns = bnet.node_sizes;
wolffd@0: onodes = find(~isemptycell(evidence));
wolffd@0: pot_type = determine_pot_type(bnet, onodes);
wolffd@0: 
wolffd@0: % prime each local kernel with evidence (if any)
wolffd@0: local_kernel = cell(1, ndoms);
wolffd@0: for i=1:ndoms
wolffd@0:   if myismember(i, exclude)
wolffd@0:     local_kernel{i} =  mk_initial_pot(pot_type, gdl.doms{i}, ns, bnet.cnodes, onodes);
wolffd@0:   else
wolffd@0:     e = bnet.equiv_class(i);
wolffd@0:     local_kernel{i} =  convert_to_pot(bnet.CPD{e}, pot_type, gdl.doms{i}(:), evidence);
wolffd@0:   end
wolffd@0: end
wolffd@0:   
wolffd@0: % initialise all msgs to 1s
wolffd@0: msg = cell(ndoms, ndoms);
wolffd@0: for i=1:ndoms
wolffd@0:   nbrs = gdl.nbrs{i};
wolffd@0:   for j=nbrs(:)'
wolffd@0:     dom = gdl.sepset{i,j};
wolffd@0:     msg{i,j} = mk_initial_pot(pot_type, dom, ns, bnet.cnodes, onodes);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: switch engine.protocol
wolffd@0:  case 'parallel', 
wolffd@0:    [engine.marginal_domains, niter] = parallel_protocol(engine, evidence, pot_type, local_kernel, msg);
wolffd@0:  case 'tree',
wolffd@0:   engine.marginal_domains = serial_protocol(engine, evidence, pot_type, local_kernel, msg);
wolffd@0:   niter = 1;
wolffd@0: end
wolffd@0: engine.niter = niter;
wolffd@0: 
wolffd@0: %fprintf('just finished %d iterations of belprop\n', niter);
wolffd@0: 
wolffd@0: if ~isempty(engine.filename)
wolffd@0:   fclose(engine.fid);
wolffd@0: end