diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine')
21 files changed, 1546 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Entries new file mode 100644 index 00000000..84a6daa1 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Entries @@ -0,0 +1,6 @@ +/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/enter_soft_ev.m/1.1.1.1/Wed May 29 15:59:56 2002// +/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002// +/pearl_dbn_inf_engine.m/1.1.1.1/Wed May 29 15:59:56 2002// +D/Old//// +D/private//// diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Repository new file mode 100644 index 00000000..b7a44128 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/dynamic/@pearl_dbn_inf_engine diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Entries new file mode 100644 index 00000000..d729c48f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Entries @@ -0,0 +1,8 @@ +/correct_smooth.m/1.1.1.1/Wed May 29 15:59:56 2002// +/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/filter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/filter_evidence_obj_oriented.m/1.1.1.1/Wed May 29 15:59:56 2002// +/smooth_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/smooth_evidence_fast.m/1.1.1.1/Wed May 29 15:59:56 2002// +/wrong_smooth.m/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Repository new file mode 100644 index 00000000..db9771c6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/correct_smooth.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/correct_smooth.m new file mode 100644 index 00000000..275afd41 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/correct_smooth.m @@ -0,0 +1,244 @@ +function [marginal, msg, loglik] = smooth_evidence(engine, evidence) +% [marginal, msg, loglik] = smooth_evidence(engine, evidence) (pearl_dbn) + +disp('warning: broken'); + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +bnet2 = dbn_to_bnet(bnet, T); +ns = bnet2.node_sizes; +hnodes = mysetdiff(1:ss, engine.onodes); +hnodes = hnodes(:)'; + +onodes2 = unroll_set(engine.onodes(:), ss, T); +onodes2 = onodes2(:)'; + +hnodes2 = unroll_set(hnodes(:), ss, T); +hnodes2 = hnodes2(:)'; + +[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2); + +msg = init_msgs(bnet2.dag, ns, evidence, bnet2.equiv_class, bnet2.CPD); + +verbose = 0; + +niter = 1; +for iter=1:niter + % FORWARD + for t=1:T + if verbose, fprintf('t=%d\n', t); end + % observed leaves send lambda to parents + for i=engine.onodes(:)' + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + for p=ps(:)' + j = engine.child_index{p}(n); % n is p's j'th child + if t > 1 + e = bnet.equiv_class(i, 2); + else + e = bnet.equiv_class(i, 1); + end + lam_msg = normalise(compute_lambda_msg(bnet.CPD{e}, n, ps, msg, p)); + msg{p}.lambda_from_child{j} = lam_msg; + if verbose, fprintf('%d sends lambda to %d\n', n, p); disp(lam_msg); end + end + end + + % update pi + for i=hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + if t==1 + e = bnet.equiv_class(i,1); + else + e = bnet.equiv_class(i,2); + end + msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg); + if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end + end + + % send pi msg to children + for i=hnodes + n = i + (t-1)*ss; + %cs = myintersect(children(bnet2.dag, n), hnodes2); + cs = children(bnet2.dag, n); + for c=cs(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns)); + msg{c}.pi_from_parent{j} = pi_msg; + if verbose, fprintf('%d sends pi to %d\n', n, c); disp(pi_msg); end + end + end + end + + % BACKWARD + for t=T:-1:1 + if verbose, fprintf('t = %d\n', t); end + % update lambda + for i=hnodes + n = i + (t-1)*ss; + cs = children(bnet2.dag, n); + msg{n}.lambda = compute_lambda(n, cs, msg, ns); + if verbose, fprintf('%d computes lambda\n', n); disp(msg{n}.lambda); end + end + % send lambda msgs to parents + for i=hnodes + n = i + (t-1)*ss; + %ps = myintersect(parents(bnet2.dag, n), hnodes2); + ps = parents(bnet2.dag, n); + for p=ps(:)' + j = engine.child_index{p}(n); % n is p's j'th child + if t > 1 + e = bnet.equiv_class(i, 2); + else + e = bnet.equiv_class(i, 1); + end + lam_msg = normalise(compute_lambda_msg(bnet.CPD{e}, n, ps, msg, p)); + msg{p}.lambda_from_child{j} = lam_msg; + if verbose, fprintf('%d sends lambda to %d\n', n, p); disp(lam_msg); end + end + end + end + +end + + +marginal = cell(ss,T); +lik = zeros(1,ss*T); +for t=1:T + for i=1:ss + n = i + (t-1)*ss; + [bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda); + marginal{i,t} = bel; + end +end + +loglik = sum(log(lik)); + + + +%%%%%%% + +function lambda = compute_lambda(n, cs, msg, ns) +% Pearl p183 eq 4.50 +lambda = prod_lambda_msgs(n, cs, msg, ns); + +%%%%%%% + +function pi_msg = compute_pi_msg(n, cs, msg, c, ns) +% Pearl p183 eq 4.53 and 4.51 +pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c); + +%%%%%%%%% + +function lam = prod_lambda_msgs(n, cs, msg, ns, except) + +if nargin < 5, except = -1; end + +lam = msg{n}.lambda_from_self(:); +lam = ones(ns(n), 1); +for i=1:length(cs) + c = cs(i); + if c ~= except + lam = lam .* msg{n}.lambda_from_child{i}; + end +end + + +%%%%%%%%% + +function msg = init_msgs(dag, ns, evidence, eclass, CPD) +% INIT_MSGS Initialize the lambda/pi message and state vectors (pearl_dbn) +% msg = init_msgs(dag, ns, evidence) + +N = length(dag); +msg = cell(1,N); +observed = ~isemptycell(evidence(:)); + +for n=1:N + ps = parents(dag, n); + msg{n}.pi_from_parent = cell(1, length(ps)); + for i=1:length(ps) + p = ps(i); + msg{n}.pi_from_parent{i} = ones(ns(p), 1); + end + + cs = children(dag, n); + msg{n}.lambda_from_child = cell(1, length(cs)); + for i=1:length(cs) + c = cs(i); + msg{n}.lambda_from_child{i} = ones(ns(n), 1); + end + + msg{n}.lambda = ones(ns(n), 1); + msg{n}.lambda_from_self = ones(ns(n), 1); + msg{n}.pi = ones(ns(n), 1); + + % Initialize the lambdas with any evidence + if observed(n) + v = evidence{n}; + %msg{n}.lambda_from_self = zeros(ns(n), 1); + %msg{n}.lambda_from_self(v) = 1; % delta function + msg{n}.lambda = zeros(ns(n), 1); + msg{n}.lambda(v) = 1; % delta function + end + +end + + +%%%%%%%% + +function msg = init_ev_msgs(engine, evidence, msg) + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +pot_type = 'd'; +t = 1; +hnodes = mysetdiff(1:ss, engine.onodes); +for i=engine.onodes(:)' + fam = family(bnet.dag, i); + e = bnet.equiv_class(i, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + temp = pot_to_marginal(CPDpot); + msg{i}.lambda_from_self = temp.T; +end +for t=2:T + for i=engine.onodes(:)' + fam = family(bnet.dag, i, 2); % extract from slice t + e = bnet.equiv_class(i, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + temp = pot_to_marginal(CPDpot); + n = i + (t-1)*ss; + msg{n}.lambda_from_self = temp.T; + end +end + + +%%%%%%%%%%% + +function msg = init_ev_msgs2(engine, evidence, msg) + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +pot_type = 'd'; +t = 1; +hnodes = mysetdiff(1:ss, engine.onodes); +for i=engine.onodes(:)' + fam = family(bnet.dag, i); + e = bnet.equiv_class(i, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + temp = pot_to_marginal(CPDpot); + msg{i}.lambda_from_self = temp.T; +end +for t=2:T + for i=engine.onodes(:)' + fam = family(bnet.dag, i, 2); % extract from slice t + e = bnet.equiv_class(i, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + temp = pot_to_marginal(CPDpot); + n = i + (t-1)*ss; + msg{n}.lambda_from_self = temp.T; + end +end + + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/enter_evidence.m new file mode 100644 index 00000000..18e7519b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/enter_evidence.m @@ -0,0 +1,123 @@ +function [engine, loglik] = enter_evidence(engine, evidence, filter) +% ENTER_EVIDENCE Add the specified evidence to the network (pearl_dbn) +% [engine, loglik] = enter_evidence(engine, evidence, filter) +% +% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector) +% If filter = 1, we do filtering, otherwise smoothing (default). + +if nargin < 3, filter = 0; end + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +bnet2 = dbn_to_bnet(bnet, T); +ns = bnet2.node_sizes; +hnodes = mysetdiff(1:ss, engine.onodes); +hnodes = hnodes(:)'; + +[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2); + +msg = init_msgs(bnet2.dag, ns, evidence); +msg = init_ev_msgs(engine, evidence, msg); + +niter = 1; +for iter=1:niter + % FORWARD + for t=1:T + % update pi + for i=1:ss %hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + if t==1 + e = bnet.equiv_class(i,1); + else + e = bnet.equiv_class(i,2); + end + msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg); + %msg{n}.pi = normalise(msg{n}.pi(:) .* msg{n}.lambda_from_self(:)); + end + % send pi msg to children + for i=1:ss % hnodes + n = i + (t-1)*ss; + cs = children(bnet2.dag, n); + for c=cs(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + msg{c}.pi_from_parent{j} = normalise(compute_pi_msg(n, cs, msg, c, ns)); + end + end + end + + if filter + disp('skipping smoothing'); + break; + end + + % BACKWARD + for t=T:-1:1 + % update lambda + for i=1:ss % hnodes + n = i + (t-1)*ss; + cs = children(bnet2.dag, n); + msg{n}.lambda = compute_lambda(n, cs, msg, ns); + end + % send lambda msgs to parents + for i=1:ss % hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + for p=ps(:)' + j = engine.child_index{p}(n); % n is p's j'th child + if t > 1 + e = bnet.equiv_class(i, 2); + else + e = bnet.equiv_class(i, 1); + end + msg{p}.lambda_from_child{j} = normalise(compute_lambda_msg(bnet.CPD{e}, n, ps, msg, p)); + end + end + end + +end + + +engine.marginal = cell(ss,T); +lik = zeros(1,ss*T); +for t=1:T + for i=1:ss + n = i + (t-1)*ss; + [bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda); + engine.marginal{i,t} = bel; + end +end + +engine.evidence = evidence; % needed by marginal_nodes and marginal_family +engine.msg = msg; % needed by marginal_family +loglik = sum(log(lik)); + + + +%%%%%%% + +function lambda = compute_lambda(n, cs, msg, ns) +% Pearl p183 eq 4.50 +lambda = prod_lambda_msgs(n, cs, msg, ns); + +%%%%%%% + +function pi_msg = compute_pi_msg(n, cs, msg, c, ns) +% Pearl p183 eq 4.53 and 4.51 +pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c); + +%%%%%%%%% + +function lam = prod_lambda_msgs(n, cs, msg, ns, except) + +if nargin < 5, except = -1; end + +lam = msg{n}.lambda_from_self(:); +%lam = ones(ns(n), 1); +for i=1:length(cs) + c = cs(i); + if c ~= except + lam = lam .* msg{n}.lambda_from_child{i}; + end +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence.m new file mode 100644 index 00000000..a3462437 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence.m @@ -0,0 +1,146 @@ +function [marginal, msg, loglik] = filter_evidence(engine, evidence) + +error('broken'); + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +onodes = engine.onodes; +hnodes = mysetdiff(1:ss, onodes); +hnodes = hnodes(:)'; + +ns = bnet.node_sizes(:); +onodes2 = [onodes(:); onodes(:)+ss]; +ns(onodes2) = 1; + +verbose = 1; +if verbose, fprintf('\nnew filtering\n'); end + +pot_type = 'd'; +niter = engine.max_iter; + +% msg(i1,t1,i2,j2) (i1,t1) -> (i2,t2) +%lambda_msg = cell(ss,T,ss,T); +%pi_msg = cell(ss,T,ss,T); + +% intra_lambda_msg(i,j,t) (i,t) -> (j,t), i is child +% inter_lambda_msg(i,j,t) (i,t+1) -> (j,t), i is child +% inter_pi_msg(i,j,t) (i,t-1) -> (j,t), i is parent +intra_lambda_msg = cell(ss,ss,T); +inter_lambda_msg = cell(ss,ss,T); +inter_pi_msg = cell(ss,ss,T); + +lambda = cell(ss,T); +pi = cell(ss,T); + +for t=1:T + for i=1:ss + lambda{i,t} = ones(ns(i), 1); + pi{i,t} = ones(ns(i), 1); + + cs = children(bnet.intra, i); + for c=cs(:)' + intra_lambda_msg{c,i,t} = ones(ns(i),1); + end + + cs = children(bnet.inter, i); + for c=cs(:)' + inter_lambda_msg{c,i,t} = ones(ns(i),1); + end + + ps = parents(bnet.inter, i); + for p=ps(:)' + inter_pi_msg{p,i,t} = ones(ns(i), 1); % not used for t==1 + end + end +end + +% each hidden node absorbs lambda from its observed child (if any) +for t=1:T + for i=hnodes + c = engine.obschild(i); + if c > 0 + if t==1 + fam = family(bnet.dag, c); + e = bnet.equiv_class(c, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + else + fam = family(bnet.dag, c, 2); % within 2 slice network + e = bnet.equiv_class(c, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + end + temp = pot_to_marginal(CPDpot); + lam_msg = normalise(temp.T); + %if verbose, fprintf('(%d,%d) sends lambda to (%d,%d)\n', c,t, i,t); disp(lam_msg); end + intra_lambda_msg{c,i,t} = lam_msg; + end + end +end + +% FORWARD +for t=1:T + % update pi + for i=hnodes + if t==1 + e = bnet.equiv_class(i,1); + temp = struct(bnet.CPD{e}); + pi{i,t} = temp.CPT; + else + e = bnet.equiv_class(i,2); + temp = struct(bnet.CPD{e}); + ps = parents(bnet.inter, i); + dom = [ps i+ss]; + pot = dpot(dom, ns(dom), temp.CPT); + for p=ps(:)' + temp = dpot(p, ns(p), inter_pi_msg{p,i,t}); + pot = multiply_by_pot(pot, temp); + end + pot = marginalize_pot(pot, i+ss); + temp = pot_to_marginal(pot); + pi{i,t} = temp.T; + %if verbose, fprintf('(%d,%d) computes pi\n', i,t); disp(pi{i,t}); end + end + + c = engine.obschild(i); + if c > 0 + pi{i,t} = normalise(pi{i,t} .* intra_lambda_msg{c,i,t}); + end + %if verbose, fprintf('(%d,%d) recomputes pi\n', i,t); disp(pi{i,t}); end + if verbose, fprintf('%d recomputes pi\n', i+(t-1)*ss); disp(pi{i,t}); end + end + + % send pi msg to children + for i=hnodes + cs = children(bnet.inter, i); + for c=cs(:)' + pot = pi{i,t}; + for k=cs(:)' + if k ~= c + pot = pot .* inter_lambda_msg{k,i,t}; + end + end + cs2 = children(bnet.intra, i); + for k=cs2(:)' + pot = pot .* intra_lambda_msg{k,i,t}; + end + pot = normalise(pot); + %if verbose, fprintf('(%d,%d) sends pi to (%d,%d)\n', i,t, c,t+1); disp(pot); end + if verbose, fprintf('%d sends pi to %d\n', i+(t-1)*ss, c+t*ss); disp(pot); end + inter_pi_msg{i,c,t+1} = pot; + end + end +end + + +marginal = cell(ss,T); +for t=1:T + for i=hnodes + %marginal{i,t} = normalise(pi{i,t} .* lambda{i,t}); + marginal{i,t} = normalise(pi{i,t}); + end +end + +loglik = 0; + +msg.inter_pi_msg = inter_pi_msg; +msg.inter_lambda_msg = inter_lambda_msg; +msg.intra_lambda_msg = intra_lambda_msg; diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence_obj_oriented.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence_obj_oriented.m new file mode 100644 index 00000000..fec80b11 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence_obj_oriented.m @@ -0,0 +1,158 @@ +function [marginal, msg, loglik] = filter_evidence_old(engine, evidence) +% [marginal, msg, loglik] = filter_evidence(engine, evidence) (pearl_dbn) + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +bnet2 = dbn_to_bnet(bnet, T); +ns = bnet2.node_sizes; +hnodes = mysetdiff(1:ss, engine.onodes); +hnodes = hnodes(:)'; + +[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2); + +msg = init_msgs(bnet2.dag, ns, evidence); +msg = init_ev_msgs(engine, evidence, msg); + +verbose = 1; +if verbose, fprintf('\nold filtering\n'); end + +for t=1:T + % update pi + for i=hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + if t==1 + e = bnet.equiv_class(i,1); + else + e = bnet.equiv_class(i,2); + end + msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg); + %if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end + msg{n}.pi = normalise(msg{n}.pi(:) .* msg{n}.lambda_from_self(:)); + if verbose, fprintf('%d recomputes pi\n', n); disp(msg{n}.pi); end + end + % send pi msg to children + for i=hnodes + n = i + (t-1)*ss; + cs = children(bnet2.dag, n); + for c=cs(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns)); + msg{c}.pi_from_parent{j} = pi_msg; + if verbose, fprintf('%d sends pi to %d\n', n,c); disp(pi_msg); end + end + end +end + + +marginal = cell(ss,T); +lik = zeros(1,ss*T); +for t=1:T + for i=1:ss + n = i + (t-1)*ss; + %[bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda); + [bel, lik(n)] = normalise(msg{n}.pi); + marginal{i,t} = bel; + end +end + +loglik = sum(log(lik)); + + + +%%%%%%% + +function lambda = compute_lambda(n, cs, msg, ns) +% Pearl p183 eq 4.50 +lambda = prod_lambda_msgs(n, cs, msg, ns); + +%%%%%%% + +function pi_msg = compute_pi_msg(n, cs, msg, c, ns) +% Pearl p183 eq 4.53 and 4.51 +pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c); + +%%%%%%%%% + +function lam = prod_lambda_msgs(n, cs, msg, ns, except) + +if nargin < 5, except = -1; end + +%lam = msg{n}.lambda_from_self(:); +lam = ones(ns(n), 1); +for i=1:length(cs) + c = cs(i); + if c ~= except + lam = lam .* msg{n}.lambda_from_child{i}; + end +end + + +%%%%%%%%%%% + +function msg = init_msgs(dag, ns, evidence) +% INIT_MSGS Initialize the lambda/pi message and state vectors (pearl_dbn) +% msg = init_msgs(dag, ns, evidence) +% +% We assume all the hidden nodes are discrete. + +N = length(dag); +msg = cell(1,N); +observed = ~isemptycell(evidence(:)); + +for n=1:N + ps = parents(dag, n); + msg{n}.pi_from_parent = cell(1, length(ps)); + for i=1:length(ps) + p = ps(i); + msg{n}.pi_from_parent{i} = ones(ns(p), 1); + end + + cs = children(dag, n); + msg{n}.lambda_from_child = cell(1, length(cs)); + for i=1:length(cs) + c = cs(i); + msg{n}.lambda_from_child{i} = ones(ns(n), 1); + end + + msg{n}.lambda = ones(ns(n), 1); + msg{n}.pi = ones(ns(n), 1); + + msg{n}.lambda_from_self = ones(ns(n), 1); +end + + +%%%%%%%%% + +function msg = init_ev_msgs(engine, evidence, msg) +% Initialize the lambdas with any evidence + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +pot_type = 'd'; +t = 1; +hnodes = mysetdiff(1:ss, engine.onodes); +for i=hnodes(:)' + c = engine.obschild(i); + if c > 0 + fam = family(bnet.dag, c); + e = bnet.equiv_class(c, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + temp = pot_to_marginal(CPDpot); + n = i; + msg{n}.lambda_from_self = temp.T; + end +end +for t=2:T + for i=hnodes(:)' + c = engine.obschild(i); + if c > 0 + fam = family(bnet.dag, c, 2); + e = bnet.equiv_class(c, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + temp = pot_to_marginal(CPDpot); + n = i + (t-1)*ss; + msg{n}.lambda_from_self = temp.T; + end + end +end diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/smooth_evidence.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/smooth_evidence.m new file mode 100644 index 00000000..554b579f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/smooth_evidence.m @@ -0,0 +1,181 @@ +function [marginal, msg, loglik] = smooth_evidence(engine, evidence) + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +onodes = engine.onodes; +hnodes = mysetdiff(1:ss, onodes); +hnodes = hnodes(:)'; + +ns = bnet.node_sizes(:); +onodes2 = [onodes(:); onodes(:)+ss]; +ns(onodes2) = 1; + +verbose = 0; +pot_type = 'd'; +niter = engine.max_iter; + +if verbose, fprintf('new smooth\n'); end + +% msg(i1,t1,i2,j2) (i1,t1) -> (i2,t2) +%lambda_msg = cell(ss,T,ss,T); +%pi_msg = cell(ss,T,ss,T); + +% intra_lambda_msg(i,j,t) (i,t) -> (j,t), i is child +% inter_lambda_msg(i,j,t) (i,t+1) -> (j,t), i is child +% inter_pi_msg(i,j,t) (i,t-1) -> (j,t), i is parent +intra_lambda_msg = cell(ss,ss,T); +inter_lambda_msg = cell(ss,ss,T); +inter_pi_msg = cell(ss,ss,T); + +lambda = cell(ss,T); +pi = cell(ss,T); + +for t=1:T + for i=1:ss + lambda{i,t} = ones(ns(i), 1); + pi{i,t} = ones(ns(i), 1); + + cs = children(bnet.intra, i); + for c=cs(:)' + intra_lambda_msg{c,i,t} = ones(ns(i),1); + end + + cs = children(bnet.inter, i); + for c=cs(:)' + inter_lambda_msg{c,i,t} = ones(ns(i),1); + end + + ps = parents(bnet.inter, i); + for p=ps(:)' + inter_pi_msg{p,i,t} = ones(ns(i), 1); % not used for t==1 + end + end +end + + +% each hidden node absorbs lambda from its observed child (if any) +for t=1:T + for i=hnodes + c = engine.obschild(i); + if c > 0 + if t==1 + fam = family(bnet.dag, c); + e = bnet.equiv_class(c, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + else + fam = family(bnet.dag, c, 2); % within 2 slice network + e = bnet.equiv_class(c, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + end + temp = pot_to_marginal(CPDpot); + lam_msg = normalise(temp.T); + intra_lambda_msg{c,i,t} = lam_msg; + end + end +end + +for iter=1:engine.max_iter + % FORWARD + for t=1:T + % update pi + for i=hnodes + if t==1 + e = bnet.equiv_class(i,1); + CPD = struct(bnet.CPD{e}); + pi{i,t} = CPD.CPT; + else + e = bnet.equiv_class(i,2); + CPD = struct(bnet.CPD{e}); + ps = parents(bnet.inter, i); + dom = [ps i+ss]; + pot = dpot(dom, ns(dom), CPD.CPT); + for p=ps(:)' + temp = dpot(p, ns(p), inter_pi_msg{p,i,t}); + pot = multiply_by_pot(pot, temp); + end + pot = marginalize_pot(pot, i+ss); + temp = pot_to_marginal(pot); + pi{i,t} = temp.T; + end + if verbose, fprintf('%d updates pi\n', i+(t-1)*ss); disp(pi{i,t}); end + end + + % send pi msg to children + for i=hnodes + cs = children(bnet.inter, i); + for c=cs(:)' + pot = pi{i,t}; + for k=cs(:)' + if k ~= c + pot = pot .* inter_lambda_msg{k,i,t}; + end + end + cs2 = children(bnet.intra, i); + for k=cs2(:)' + pot = pot .* intra_lambda_msg{k,i,t}; + end + inter_pi_msg{i,c,t+1} = normalise(pot); + if verbose, fprintf('%d sends pi to %d\n', i+(t-1)*ss, c+t*ss); disp(inter_pi_msg{i,c,t+1}); end + end + end + end + + if verbose, fprintf('backwards\n'); end + % BACKWARD + for t=T:-1:1 + % update lambda + for i=hnodes + pot = ones(ns(i), 1); + cs = children(bnet.inter, i); + for c=cs(:)' + pot = pot .* inter_lambda_msg{c,i,t}; + end + cs = children(bnet.intra, i); + for c=cs(:)' + pot = pot .* intra_lambda_msg{c,i,t}; + end + lambda{i,t} = normalise(pot); + if verbose, fprintf('%d computes lambda\n', i+(t-1)*ss); disp(lambda{i,t}); end + end + + % send lambda msgs to hidden parents in prev slcie + for i=hnodes + ps = parents(bnet.inter, i); + if t > 1 + e = bnet.equiv_class(i, 2); + CPD = struct(bnet.CPD{e}); + fam = [ps i+ss]; + for p=ps(:)' + pot = dpot(fam, ns(fam), CPD.CPT); + temp = dpot(i+ss, ns(i), lambda{i,t}); + pot = multiply_by_pot(pot, temp); + for k=ps(:)' + if k ~= p + temp = dpot(k, ns(k), inter_pi_msg{k,i,t}); + pot = multiply_by_pot(pot, temp); + end + end + pot = marginalize_pot(pot, p); + temp = pot_to_marginal(pot); + inter_lambda_msg{i,p,t-1} = normalise(temp.T); + if verbose, fprintf('%d sends lambda to %d\n', i+(t-1)*ss, p+(t-2)*ss); disp(inter_lambda_msg{i,p,t-1}); end + end + end + end + end +end + + + +marginal = cell(ss,T); +for t=1:T + for i=hnodes + marginal{i,t} = normalise(pi{i,t} .* lambda{i,t}); + end +end + +loglik = 0; + +msg.inter_pi_msg = inter_pi_msg; +msg.inter_lambda_msg = inter_lambda_msg; +msg.intra_lambda_msg = intra_lambda_msg; diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/smooth_evidence_fast.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/smooth_evidence_fast.m new file mode 100644 index 00000000..8f4ebd2f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/smooth_evidence_fast.m @@ -0,0 +1,179 @@ +function [marginal, msg, loglik] = smooth_evidence_fast(engine, evidence) + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +onodes = engine.onodes; +hnodes = mysetdiff(1:ss, onodes); +hnodes = hnodes(:)'; + +ns = bnet.node_sizes(:); +onodes2 = [onodes(:); onodes(:)+ss]; +ns(onodes2) = 1; + +verbose = 0; +pot_type = 'd'; +niter = engine.max_iter; + +if verbose, fprintf('new smooth\n'); end + +% msg(i1,t1,i2,j2) (i1,t1) -> (i2,t2) +%lambda_msg = cell(ss,T,ss,T); +%pi_msg = cell(ss,T,ss,T); + +% intra_lambda_msg(i,j,t) (i,t) -> (j,t), i is child +% inter_lambda_msg(i,j,t) (i,t+1) -> (j,t), i is child +% inter_pi_msg(i,j,t) (i,t-1) -> (j,t), i is parent +intra_lambda_msg = cell(ss,ss,T); +inter_lambda_msg = cell(ss,ss,T); +inter_pi_msg = cell(ss,ss,T); + +lambda = cell(ss,T); +pi = cell(ss,T); + +for t=1:T + for i=1:ss + lambda{i,t} = ones(ns(i), 1); + pi{i,t} = ones(ns(i), 1); + + cs = children(bnet.intra, i); + for c=cs(:)' + intra_lambda_msg{c,i,t} = ones(ns(i),1); + end + + cs = children(bnet.inter, i); + for c=cs(:)' + inter_lambda_msg{c,i,t} = ones(ns(i),1); + end + + ps = parents(bnet.inter, i); + for p=ps(:)' + inter_pi_msg{p,i,t} = ones(ns(i), 1); % not used for t==1 + end + end +end + + +% each hidden node absorbs lambda from its observed child (if any) +for t=1:T + for i=hnodes + c = engine.obschild(i); + if c > 0 + if t==1 + fam = family(bnet.dag, c); + e = bnet.equiv_class(c, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + else + fam = family(bnet.dag, c, 2); % within 2 slice network + e = bnet.equiv_class(c, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + end + temp = pot_to_marginal(CPDpot); + lam_msg = normalise(temp.T); + intra_lambda_msg{c,i,t} = lam_msg; + end + end +end + +for iter=1:engine.max_iter + % FORWARD + for t=1:T + % update pi + for i=hnodes + if t==1 + e = bnet.equiv_class(i,1); + temp = struct(bnet.CPD{e}); + pi{i,t} = temp.CPT; + else + e = bnet.equiv_class(i,2); + CPD = struct(bnet.CPD{e}); + ps = parents(bnet.inter, i); + temp = CPD.CPT; + for p=ps(:)' + temp(:) = temp(:) .* inter_pi_msg{p,i,t}(engine.mult_parent_ndx{i,p}); + end + dom = [ps i+ss]; + pot = dpot(dom, ns(dom), temp); + pot = marginalize_pot(pot, i+ss); + temp = pot_to_marginal(pot); + pi{i,t} = temp.T; + end + if verbose, fprintf('%d updates pi\n', i+(t-1)*ss); disp(pi{i,t}); end + end + + % send pi msg to children + for i=hnodes + cs = children(bnet.inter, i); + for c=cs(:)' + pot = pi{i,t}; + for k=cs(:)' + if k ~= c + pot = pot .* inter_lambda_msg{k,i,t}; + end + end + cs2 = children(bnet.intra, i); + for k=cs2(:)' + pot = pot .* intra_lambda_msg{k,i,t}; + end + inter_pi_msg{i,c,t+1} = normalise(pot); + if verbose, fprintf('%d sends pi to %d\n', i+(t-1)*ss, c+t*ss); disp(inter_pi_msg{i,c,t+1}); end + end + end + end + + if verbose, fprintf('backwards\n'); end + % BACKWARD + for t=T:-1:1 + % update lambda + for i=hnodes + pot = ones(ns(i), 1); + cs = children(bnet.inter, i); + for c=cs(:)' + pot = pot .* inter_lambda_msg{c,i,t}; + end + cs = children(bnet.intra, i); + for c=cs(:)' + pot = pot .* intra_lambda_msg{c,i,t}; + end + lambda{i,t} = normalise(pot); + if verbose, fprintf('%d computes lambda\n', i+(t-1)*ss); disp(lambda{i,t}); end + end + + % send lambda msgs to hidden parents in prev slcie + for i=hnodes + ps = parents(bnet.inter, i); + if t > 1 + e = bnet.equiv_class(i, 2); + CPD = struct(bnet.CPD{e}); + for p=ps(:)' + temp = CPD.CPT(:) .* lambda{i,t}(engine.mult_self_ndx{i}); + for k=ps(:)' + if k ~= p + temp(:) = temp(:) .* inter_pi_msg{k,i,t}(engine.mult_parent_ndx{i,k}); + end + end + fam = [ps i+ss]; + pot = dpot(fam, ns(fam), temp); + pot = marginalize_pot(pot, p); + temp = pot_to_marginal(pot); + inter_lambda_msg{i,p,t-1} = normalise(temp.T); + if verbose, fprintf('%d sends lambda to %d\n', i+(t-1)*ss, p+(t-2)*ss); disp(inter_lambda_msg{i,p,t-1}); end + end + end + end + end +end + + + +marginal = cell(ss,T); +for t=1:T + for i=hnodes + marginal{i,t} = normalise(pi{i,t} .* lambda{i,t}); + end +end + +loglik = 0; + +msg.inter_pi_msg = inter_pi_msg; +msg.inter_lambda_msg = inter_lambda_msg; +msg.intra_lambda_msg = intra_lambda_msg; diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/wrong_smooth.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/wrong_smooth.m new file mode 100644 index 00000000..d66d61ad --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/Old/wrong_smooth.m @@ -0,0 +1,210 @@ +function [marginal, msg, loglik] = smooth_evidence(engine, evidence) +% [marginal, msg, loglik] = smooth_evidence(engine, evidence) (pearl_dbn) + +disp('warning: pearl_dbn smoothing is broken'); + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +bnet2 = dbn_to_bnet(bnet, T); +ns = bnet2.node_sizes; +hnodes = mysetdiff(1:ss, engine.onodes); +hnodes = hnodes(:)'; + +onodes2 = unroll_set(engine.onodes(:), ss, T); +onodes2 = onodes2(:)'; + +hnodes2 = unroll_set(hnodes(:), ss, T); +hnodes2 = hnodes2(:)'; + +[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2); + +msg = init_msgs(bnet2.dag, ns, evidence, bnet2.equiv_class, bnet2.CPD); + +verbose = 0; +pot_type = 'd'; +niter = 1; +for iter=1:niter + % FORWARD + for t=1:T + if verbose, fprintf('t=%d\n', t); end + + % each hidden node absorbs lambda from its observed child (if any) + for i=hnodes + c = engine.obschild(i); + if c > 0 + if t==1 + fam = family(bnet.dag, c); + e = bnet.equiv_class(c, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + else + fam = family(bnet.dag, 2); % within 2 slice network + e = bnet.equiv_class(c, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + end + temp = pot_to_marginal(CPDpot); + n = i + (t-1)*ss; + lam_msg = normalise(temp.T); + j = engine.child_index{n}(c+(t-1)*ss); + assert(j==1); + msg{n}.lambda_from_child{j} = lam_msg; + if verbose, fprintf('%d sends lambda to %d\n', c + (t-1)*ss, n); disp(lam_msg); end + end + end + + % update pi + for i=hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + if t==1 + e = bnet.equiv_class(i,1); + else + e = bnet.equiv_class(i,2); + end + msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg); + if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end + end + + % send pi msg to children in next slice + for i=hnodes + n = i + (t-1)*ss; + %cs = myintersect(children(bnet2.dag, n), hnodes2); + cs = children(bnet2.dag, n); + for c=cs(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns)); + msg{c}.pi_from_parent{j} = pi_msg; + if verbose, fprintf('%d sends pi to %d\n', n, c); disp(pi_msg); end + end + end + end + + % BACKWARD + for t=T:-1:1 + if verbose, fprintf('t = %d\n', t); end + + % update lambda + for i=hnodes + n = i + (t-1)*ss; + cs = children(bnet2.dag, n); + msg{n}.lambda = compute_lambda(n, cs, msg, ns); + if verbose, fprintf('%d computes lambda\n', n); disp(msg{n}.lambda); end + end + + % send lambda msgs to hidden parents in prev slcie + for i=hnodes + n = i + (t-1)*ss; + %ps = myintersect(parents(bnet2.dag, n), hnodes2); + ps = parents(bnet2.dag, n); + for p=ps(:)' + j = engine.child_index{p}(n); % n is p's j'th child + if t > 1 + e = bnet.equiv_class(i, 2); + else + e = bnet.equiv_class(i, 1); + end + lam_msg = normalise(compute_lambda_msg(bnet.CPD{e}, n, ps, msg, p)); + msg{p}.lambda_from_child{j} = lam_msg; + if verbose, fprintf('%d sends lambda to %d\n', n, p); disp(lam_msg); end + end + end + + % send pi msg to observed children + if 0 + for i=hnodes + n = i + (t-1)*ss; + cs = myintersect(children(bnet2.dag, n), onodes2); + %cs = children(bnet2.dag, n); + for c=cs(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns)); + msg{c}.pi_from_parent{j} = pi_msg; + if verbose, fprintf('%d sends pi to %d\n', n, c); disp(pi_msg); end + end + end + end + + end +end + + +marginal = cell(ss,T); +lik = zeros(1,ss*T); +for t=1:T + for i=hnodes + n = i + (t-1)*ss; + [bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda); + marginal{i,t} = bel; + end +end + +loglik = 0; +%loglik = sum(log(lik)); + + + +%%%%%%% + +function lambda = compute_lambda(n, cs, msg, ns) +% Pearl p183 eq 4.50 +lambda = prod_lambda_msgs(n, cs, msg, ns); + +%%%%%%% + +function pi_msg = compute_pi_msg(n, cs, msg, c, ns) +% Pearl p183 eq 4.53 and 4.51 +pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c); + +%%%%%%%%% + +function lam = prod_lambda_msgs(n, cs, msg, ns, except) + +if nargin < 5, except = -1; end + +%lam = msg{n}.lambda_from_self(:); +lam = ones(ns(n), 1); +for i=1:length(cs) + c = cs(i); + if c ~= except + lam = lam .* msg{n}.lambda_from_child{i}; + end +end + + +%%%%%%%%% + +function msg = init_msgs(dag, ns, evidence, eclass, CPD) +% INIT_MSGS Initialize the lambda/pi message and state vectors (pearl_dbn) +% msg = init_msgs(dag, ns, evidence) + +N = length(dag); +msg = cell(1,N); +observed = ~isemptycell(evidence(:)); + +for n=1:N + ps = parents(dag, n); + msg{n}.pi_from_parent = cell(1, length(ps)); + for i=1:length(ps) + p = ps(i); + msg{n}.pi_from_parent{i} = ones(ns(p), 1); + end + + cs = children(dag, n); + msg{n}.lambda_from_child = cell(1, length(cs)); + for i=1:length(cs) + c = cs(i); + msg{n}.lambda_from_child{i} = ones(ns(n), 1); + end + + msg{n}.lambda = ones(ns(n), 1); + msg{n}.pi = ones(ns(n), 1); + + % Initialize the lambdas with any evidence + if observed(n) + v = evidence{n}; + msg{n}.lambda = zeros(ns(n), 1); + msg{n}.lambda(v) = 1; % delta function + msg{n}.lambda = []; + end + +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/enter_evidence.m new file mode 100644 index 00000000..624cac96 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/enter_evidence.m @@ -0,0 +1,35 @@ +function [engine, loglik] = enter_evidence(engine, evidence, varargin) +% ENTER_EVIDENCE Add the specified evidence to the network (loopy_dbn) +% [engine, loglik] = enter_evidence(engine, evidence, ....) +% +% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector) +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% maximize - if 1, does max-product (not yet supported), else sum-product [0] +% +% e.g., engine = enter_evidence(engine, ev, 'maximize', 1) + +maximize = 0; +filter = 0; + +% parse optional params +args = varargin; +nargs = length(args); +if nargs > 0 + for i=1:2:nargs + switch args{i}, + case 'maximize', maximize = args{i+1}; + case 'filter', filter = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + +assert(~maximize); +assert(~filter); + +[engine.marginal, engine.msg, loglik] = enter_soft_ev(engine, evidence); +engine.evidence = evidence; % needed by marginal_nodes and marginal_family diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/enter_soft_ev.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/enter_soft_ev.m new file mode 100644 index 00000000..5c88f53f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/enter_soft_ev.m @@ -0,0 +1,137 @@ +function [marginal, msg, loglik] = enter_soft_ev(engine, evidence) +% [marginal, msg, loglik] = smooth_evidence(engine, evidence) (pearl_dbn) + + +[ss T] = size(evidence); +bnet = bnet_from_engine(engine); +bnet2 = dbn_to_bnet(bnet, T); +ns = bnet2.node_sizes; +hnodes = mysetdiff(1:ss, engine.onodes); +hnodes = hnodes(:)'; + +onodes2 = unroll_set(engine.onodes(:), ss, T); +onodes2 = onodes2(:)'; + +hnodes2 = unroll_set(hnodes(:), ss, T); +hnodes2 = hnodes2(:)'; + +[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2); + +rand_init = 0; +use_ev = 0; +msg = init_pearl_msgs(bnet2.dag, ns, evidence, rand_init, use_ev); +msg = init_pearl_dbn_ev_msgs(bnet, evidence, engine); + +verbose = 0; +pot_type = 'd'; +niter = engine.max_iter; + +if verbose, fprintf('old smooth\n'); end + +for iter=1:niter + % FORWARD + for t=1:T + if verbose, fprintf('t=%d\n', t); end + + % update pi + for i=hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + if t==1 + e = bnet.equiv_class(i,1); + else + e = bnet.equiv_class(i,2); + end + msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg); + if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end + end + + % send pi msg to children + for i=hnodes + n = i + (t-1)*ss; + %cs = myintersect(children(bnet2.dag, n), hnodes2); + cs = children(bnet2.dag, n); % must use all children to get index right + for c=cs(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns)); + msg{c}.pi_from_parent{j} = pi_msg; + if verbose, fprintf('%d sends pi to %d\n', n, c); disp(pi_msg); end + end + end + end + + % BACKWARD + for t=T:-1:1 + if verbose, fprintf('t = %d\n', t); end + + % update lambda + for i=hnodes + n = i + (t-1)*ss; + cs = children(bnet2.dag, n); + msg{n}.lambda = compute_lambda(n, cs, msg, ns); + if verbose, fprintf('%d computes lambda\n', n); disp(msg{n}.lambda); end + end + + % send lambda msgs to hidden parents in prev slcie + for i=hnodes + n = i + (t-1)*ss; + ps = parents(bnet2.dag, n); + for p=ps(:)' + j = engine.child_index{p}(n); % n is p's j'th child + if t > 1 + e = bnet.equiv_class(i, 2); + else + e = bnet.equiv_class(i, 1); + end + lam_msg = normalise(compute_lambda_msg(bnet.CPD{e}, n, ps, msg, p)); + msg{p}.lambda_from_child{j} = lam_msg; + if verbose, fprintf('%d sends lambda to %d\n', n, p); disp(lam_msg); end + end + end + + end +end + + +marginal = cell(ss,T); +lik = zeros(1,ss*T); +for t=1:T + for i=hnodes + n = i + (t-1)*ss; + [bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda); + marginal{i,t} = bel; + end +end + +loglik = 0; +%loglik = sum(log(lik)); + + + +%%%%%%% + +function lambda = compute_lambda(n, cs, msg, ns) +% Pearl p183 eq 4.50 +lambda = prod_lambda_msgs(n, cs, msg, ns); + +%%%%%%% + +function pi_msg = compute_pi_msg(n, cs, msg, c, ns) +% Pearl p183 eq 4.53 and 4.51 +pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c); + +%%%%%%%%% + +function lam = prod_lambda_msgs(n, cs, msg, ns, except) + +if nargin < 5, except = -1; end + +%lam = msg{n}.lambda_from_self(:); +lam = ones(ns(n), 1); +for i=1:length(cs) + c = cs(i); + if c ~= except + lam = lam .* msg{n}.lambda_from_child{i}; + end +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/marginal_nodes.m new file mode 100644 index 00000000..9440459a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/marginal_nodes.m @@ -0,0 +1,18 @@ +function marginal = marginal_nodes(engine, nodes, t) +% MARGINAL_NODES Compute the marginal on the specified query nodes (pearl_dbn) +% marginal = marginal_nodes(engine, i, t) +% returns Pr(X(i,t) | Y(1:T)), where X(i,t) is the i'th node in the t'th slice. +% If enter_evidence used filtering instead of smoothing, this will return Pr(X(i,t) | Y(1:t)). + +if nargin < 3, t = 1; end +assert(length(nodes)==1); +i = nodes(end); +if ~myismember(i, engine.onodes) + marginal.T = engine.marginal{i,t}; +else + marginal.T = 1; % observed +end + +% we convert the domain to the unrolled numbering system +% so that update_ess extracts the right evidence. +marginal.domain = nodes+(t-1)*engine.ss; diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/pearl_dbn_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/pearl_dbn_inf_engine.m new file mode 100644 index 00000000..2e0509fd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/pearl_dbn_inf_engine.m @@ -0,0 +1,65 @@ +function engine = pearl_dbn_inf_engine(bnet, varargin) +% LOOPY_DBN_INF_ENGINE Loopy Pearl version of forwards-backwards +% engine = loopy_dbn_inf_engine(bnet, ...) +% +% Optional arguments +% 'max_iter' - specifies the max num. forward-backward passes to perform [1] +% 'tol' - as in loopy_pearl [1e-3] +% 'momentum' - as in loopy_pearl [0] + +error('pearl_dbn does not work yet') + +max_iter = 1; +tol = 1e-3; +momentum = 0; + +if nargin >= 2 + args = varargin; + nargs = length(args); + for i=1:2:nargs + switch args{i}, + case 'max_iter', max_iter = args{i+1}; + case 'tol', tol = args{i+1}; + case 'momentum', momentum = args{i+1}; + end + end +end + + +engine.max_iter = max_iter; +engine.tol = tol; +engine.momentum = momentum; +engine.pearl_engine = []; +engine.T = []; +engine.ss = length(bnet.intra); + +engine.marginal = []; +engine.evidence = []; +engine.msg = []; +engine.parent_index = []; +engine.child_index = []; +%[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet); % need to unroll first + +ss = length(bnet.intra); +engines.ss = ss; +onodes = bnet.observed; +hnodes = mysetdiff(1:ss, onodes); +obschild = zeros(1,ss); +for i=hnodes(:)' + %ocs = myintersect(children(bnet.dag, i), onodes); + ocs = children(bnet.intra, i); + assert(length(ocs) <= 1); + if length(ocs)==1 + obschild(i) = ocs(1); + end +end +engine.obschild = obschild; + +engine.mult_self_ndx = []; +engine.mult_parent_ndx = []; +engine.marg_self_ndx = []; +engine.marg_parent_ndx = []; + + +engine = class(engine, 'loopy_dbn_inf_engine', inf_engine(bnet)); + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Entries new file mode 100644 index 00000000..e35b6662 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Entries @@ -0,0 +1,2 @@ +/init_pearl_dbn_ev_msgs.m/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Repository new file mode 100644 index 00000000..2cb0fa7a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/dynamic/@pearl_dbn_inf_engine/private diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/init_pearl_dbn_ev_msgs.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/init_pearl_dbn_ev_msgs.m new file mode 100644 index 00000000..893af2ae --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@pearl_dbn_inf_engine/private/init_pearl_dbn_ev_msgs.m @@ -0,0 +1,28 @@ +function msg = init_pearl_dbn_ev_msgs(bnet, evidence, engine) + +[ss T] = size(evidence); +pot_type = 'd'; + +% each hidden node absorbs lambda from its observed child (if any) +for t=1:T + for i=hnodes + c = engine.obschild(i); + if c > 0 + if t==1 + fam = family(bnet.dag, c); + e = bnet.equiv_class(c, 1); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1)); + else + fam = family(bnet.dag, c, 2); % within 2 slice network + e = bnet.equiv_class(c, 2); + CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t)); + end + temp = pot_to_marginal(CPDpot); + n = i + (t-1)*ss; + lam_msg = normalise(temp.T); + j = engine.child_index{n}(c+(t-1)*ss); + assert(j==1); + msg{n}.lambda_from_child{j} = lam_msg; + end + end +end |