Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/static/@pearl_inf_engine/enter_evidence.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [engine, loglik, niter] = enter_evidence(engine, evidence, varargin)
wolffd@0 2 % ENTER_EVIDENCE Add the specified evidence to the network (pearl)
wolffd@0 3 % [engine, loglik, num_iter] = enter_evidence(engine, evidence, ...)
wolffd@0 4 % evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)
wolffd@0 5 %
wolffd@0 6 % The following optional arguments can be specified in the form of name/value pa irs:
wolffd@0 7 % [default value in brackets]
wolffd@0 8 %
wolffd@0 9 % maximize - if 1, does max-product instead of sum-product [0]
wolffd@0 10 % 'filename' - msgs will be printed to this file, so you can assess convergence while it runs [engine.filename]
wolffd@0 11 %
wolffd@0 12 % e.g., engine = enter_evidence(engine, ev, 'maximize', 1)
wolffd@0 13 %
wolffd@0 14 % For discrete nodes, loglik is the negative Bethe free energy evaluated at the final beliefs.
wolffd@0 15 % For Gaussian nodes, loglik is currently always 0.
wolffd@0 16 %
wolffd@0 17 % 'num_iter' returns the number of iterations used.
wolffd@0 18
wolffd@0 19 maximize = 0;
wolffd@0 20 filename = engine.filename;
wolffd@0 21
wolffd@0 22 % parse optional params
wolffd@0 23 args = varargin;
wolffd@0 24 nargs = length(args);
wolffd@0 25 if nargs > 0
wolffd@0 26 for i=1:2:nargs
wolffd@0 27 switch args{i},
wolffd@0 28 case 'maximize', maximize = args{i+1};
wolffd@0 29 case 'filename', filename = args{i+1};
wolffd@0 30 otherwise,
wolffd@0 31 error(['invalid argument name ' args{i}]);
wolffd@0 32 end
wolffd@0 33 end
wolffd@0 34 end
wolffd@0 35
wolffd@0 36
wolffd@0 37 if maximize
wolffd@0 38 error('can''t handle max-prop yet')
wolffd@0 39 end
wolffd@0 40
wolffd@0 41 engine.maximize = maximize;
wolffd@0 42 engine.filename = filename;
wolffd@0 43 engine.bel = []; % reset if necessary
wolffd@0 44
wolffd@0 45 bnet = bnet_from_engine(engine);
wolffd@0 46 N = length(bnet.dag);
wolffd@0 47 ns = bnet.node_sizes(:);
wolffd@0 48
wolffd@0 49 observed_bitv = ~isemptycell(evidence);
wolffd@0 50 disconnected = find(engine.disconnected_nodes_bitv);
wolffd@0 51 if ~all(observed_bitv(disconnected))
wolffd@0 52 error(['The following discrete nodes must be observed: ' num2str(disconnected)])
wolffd@0 53 end
wolffd@0 54 msg = init_pearl_msgs(engine.msg_type, engine.msg_dag, ns, evidence);
wolffd@0 55
wolffd@0 56 niter = 1;
wolffd@0 57 switch engine.protocol
wolffd@0 58 case 'parallel', [msg, niter] = parallel_protocol(engine, evidence, msg);
wolffd@0 59 case 'tree', msg = tree_protocol(engine, evidence, msg);
wolffd@0 60 otherwise,
wolffd@0 61 error(['unrecognized protocol ' engine.protocol])
wolffd@0 62 end
wolffd@0 63 engine.niter = niter;
wolffd@0 64
wolffd@0 65 engine.marginal = cell(1,N);
wolffd@0 66 nodes = find(~engine.disconnected_nodes_bitv);
wolffd@0 67 for n=nodes(:)'
wolffd@0 68 engine.marginal{n} = compute_bel(engine.msg_type, msg{n}.pi, msg{n}.lambda);
wolffd@0 69 end
wolffd@0 70
wolffd@0 71 engine.evidence = evidence; % needed by marginal_nodes and marginal_family
wolffd@0 72 engine.msg = msg; % needed by marginal_family
wolffd@0 73
wolffd@0 74 if (nargout >= 2)
wolffd@0 75 if (engine.msg_type == 'd')
wolffd@0 76 loglik = bethe_free_energy(engine, evidence);
wolffd@0 77 else
wolffd@0 78 loglik = 0;
wolffd@0 79 end
wolffd@0 80 end
wolffd@0 81
wolffd@0 82
wolffd@0 83
wolffd@0 84 %%%%%%%%%%%
wolffd@0 85
wolffd@0 86 function msg = init_pearl_msgs(msg_type, dag, ns, evidence)
wolffd@0 87 % INIT_MSGS Initialize the lambda/pi message and state vectors
wolffd@0 88 % msg = init_msgs(dag, ns, evidence)
wolffd@0 89 %
wolffd@0 90
wolffd@0 91 N = length(dag);
wolffd@0 92 msg = cell(1,N);
wolffd@0 93 observed = ~isemptycell(evidence);
wolffd@0 94 lam_msg = 1;
wolffd@0 95
wolffd@0 96 for n=1:N
wolffd@0 97 ps = parents(dag, n);
wolffd@0 98 msg{n}.pi_from_parent = cell(1, length(ps));
wolffd@0 99 for i=1:length(ps)
wolffd@0 100 p = ps(i);
wolffd@0 101 msg{n}.pi_from_parent{i} = mk_msg(msg_type, ns(p));
wolffd@0 102 end
wolffd@0 103
wolffd@0 104 cs = children(dag, n);
wolffd@0 105 msg{n}.lambda_from_child = cell(1, length(cs));
wolffd@0 106 for i=1:length(cs)
wolffd@0 107 c = cs(i);
wolffd@0 108 msg{n}.lambda_from_child{i} = mk_msg(msg_type, ns(n), lam_msg);
wolffd@0 109 end
wolffd@0 110
wolffd@0 111 msg{n}.lambda = mk_msg(msg_type, ns(n), lam_msg);
wolffd@0 112 msg{n}.pi = mk_msg(msg_type, ns(n));
wolffd@0 113
wolffd@0 114 if observed(n)
wolffd@0 115 msg{n}.lambda_from_self = mk_msg_with_evidence(msg_type, ns(n), evidence{n});
wolffd@0 116 else
wolffd@0 117 msg{n}.lambda_from_self = mk_msg(msg_type, ns(n), lam_msg);
wolffd@0 118 end
wolffd@0 119 end
wolffd@0 120
wolffd@0 121
wolffd@0 122
wolffd@0 123 %%%%%%%%%
wolffd@0 124
wolffd@0 125 function msg = mk_msg(msg_type, sz, is_lambda_msg)
wolffd@0 126
wolffd@0 127 if nargin < 3, is_lambda_msg = 0; end
wolffd@0 128
wolffd@0 129 switch msg_type
wolffd@0 130 case 'd', msg = ones(sz, 1);
wolffd@0 131 case 'g',
wolffd@0 132 if is_lambda_msg
wolffd@0 133 msg.precision = zeros(sz, sz);
wolffd@0 134 msg.info_state = zeros(sz, 1);
wolffd@0 135 else
wolffd@0 136 msg.Sigma = zeros(sz, sz);
wolffd@0 137 msg.mu = zeros(sz,1);
wolffd@0 138 end
wolffd@0 139 end
wolffd@0 140
wolffd@0 141 %%%%%%%%%%%%
wolffd@0 142
wolffd@0 143 function msg = mk_msg_with_evidence(msg_type, sz, val)
wolffd@0 144
wolffd@0 145 switch msg_type
wolffd@0 146 case 'd',
wolffd@0 147 msg = zeros(sz, 1);
wolffd@0 148 msg(val) = 1;
wolffd@0 149 case 'g',
wolffd@0 150 %msg.observed_val = val(:);
wolffd@0 151 msg.precision = inf;
wolffd@0 152 msg.mu = val(:);
wolffd@0 153 end