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

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