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
|