diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/general')
64 files changed, 2654 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Entries b/sourcecodes/bnt-master/BNT/general/CVS/Entries new file mode 100644 index 00000000..de13cc9a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/CVS/Entries @@ -0,0 +1,50 @@ +/add_ev_to_dmarginal.m/1.1.1.1/Thu Jun 27 20:34:32 2002// +/add_evidence_to_gmarginal.m/1.1.1.1/Wed May 29 15:59:56 2002// +/bnet_to_fgraph.m/1.1.1.1/Wed May 29 15:59:54 2002// +/compute_fwd_interface.m/1.1.1.1/Wed May 29 15:59:54 2002// +/compute_interface_nodes.m/1.1.1.1/Wed May 29 15:59:54 2002// +/compute_joint_pot.m/1.1.1.1/Mon Jun 7 15:50:34 2004// +/compute_minimal_interface.m/1.1.1.1/Wed May 29 15:59:54 2002// +/convert_dbn_CPDs_to_pots.m/1.1.1.1/Fri Nov 22 22:35:00 2002// +/convert_dbn_CPDs_to_tables.m/1.1.1.1/Thu Jan 23 18:44:50 2003// +/convert_dbn_CPDs_to_tables1.m/1.1.1.1/Thu Jan 23 18:49:48 2003// +/convert_dbn_CPDs_to_tables_slow.m/1.1.1.1/Wed May 29 15:59:58 2002// +/dbn_to_bnet.m/1.1.1.1/Wed May 29 15:59:54 2002// +/dbn_to_hmm.m/1.1.1.1/Sun Feb 2 00:23:38 2003// +/determine_elim_constraints.m/1.1.1.1/Wed May 29 15:59:54 2002// +/dispcpt.m/1.1.1.1/Wed May 29 15:59:58 2002// +/do_intervention.m/1.1.1.1/Wed May 29 15:59:54 2002// +/dsep.m/1.1.1.1/Wed May 29 15:59:54 2002// +/dsep_test.m/1.1.1.1/Sat Jan 18 23:10:16 2003// +/enumerate_scenarios.m/1.1.1.1/Wed May 29 15:59:54 2002// +/fgraph_to_bnet.m/1.1.1.1/Wed May 29 15:59:54 2002// +/hodbn_to_bnet.m/1.1.1.1/Wed Jul 24 14:48:06 2002// +/is_mnet.m/1.1.1.1/Sun Jun 16 20:01:22 2002// +/linear_gaussian_to_cpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/log_lik_complete.m/1.1.1.1/Wed May 29 15:59:54 2002// +/log_marg_lik_complete.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_bnet.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_dbn.m/1.1.1.1/Sat Feb 1 19:42:14 2003// +/mk_fgraph.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_fgraph_given_ev.m/1.1.1.1/Mon Jun 24 18:56:26 2002// +/mk_higher_order_dbn.m/1.1.1.1/Tue Jul 23 13:17:04 2002// +/mk_limid.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_mnet.m/1.1.1.1/Sun Jun 16 19:52:12 2002// +/mk_mrf2.m/1.1.1.1/Tue Dec 31 22:06:48 2002// +/mk_mutilated_samples.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_named_CPT.m/1.1.1.1/Tue Mar 30 17:18:54 2004// +/mk_slice_and_half_dbn.m/1.1.1.1/Wed May 29 15:59:54 2002// +/noisyORtoTable.m/1.1.1.1/Mon Aug 2 22:05:58 2004// +/partition_dbn_nodes.m/1.1.1.1/Wed May 29 15:59:54 2002// +/partition_matrix_vec_3.m/1.1.1.1/Wed May 29 15:59:58 2002// +/sample_bnet.m/1.1.1.1/Thu Jun 10 01:17:46 2004// +/sample_bnet_nocell.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_dbn.m/1.1.1.1/Wed May 29 15:59:54 2002// +/score_bnet_complete.m/1.1.1.1/Wed May 29 15:59:54 2002// +/shrink_obs_dims_in_gaussian.m/1.1.1.1/Wed May 29 15:59:58 2002// +/shrink_obs_dims_in_table.m/1.1.1.1/Wed May 29 15:59:58 2002// +/solve_limid.m/1.1.1.1/Mon Jun 7 15:48:02 2004// +/unroll_dbn_topology.m/1.1.1.1/Wed May 29 15:59:54 2002// +/unroll_higher_order_topology.m/1.1.1.1/Fri May 31 10:25:58 2002// +/unroll_set.m/1.1.1.1/Mon Dec 16 17:57:14 2002// +D diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/general/CVS/Entries.Log new file mode 100644 index 00000000..24f16336 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/CVS/Entries.Log @@ -0,0 +1 @@ +A D/Old//// diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Repository b/sourcecodes/bnt-master/BNT/general/CVS/Repository new file mode 100644 index 00000000..08dc139c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/general diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Root b/sourcecodes/bnt-master/BNT/general/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/general/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/general/Old/CVS/Entries new file mode 100644 index 00000000..a1785ca4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/CVS/Entries @@ -0,0 +1,9 @@ +/bnet_to_gdl_graph.m/1.1.1.1/Wed May 29 15:59:54 2002// +/calc_mpe.m/1.1.1.1/Mon Jun 17 21:58:38 2002// +/calc_mpe_bucket.m/1.1.1.1/Wed May 29 15:59:54 2002// +/calc_mpe_dbn.m/1.1.1.1/Wed May 29 15:59:54 2002// +/calc_mpe_given_inf_engine.m/1.1.1.1/Wed May 29 15:59:54 2002// +/calc_mpe_global.m/1.1.1.1/Wed May 29 15:59:54 2002// +/compute_interface_nodes.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_gdl_graph.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/general/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/general/Old/CVS/Repository new file mode 100644 index 00000000..46960c1c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/general/Old diff --git a/sourcecodes/bnt-master/BNT/general/Old/CVS/Root b/sourcecodes/bnt-master/BNT/general/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/general/Old/bnet_to_gdl_graph.m b/sourcecodes/bnt-master/BNT/general/Old/bnet_to_gdl_graph.m new file mode 100644 index 00000000..d6ff45f3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/bnet_to_gdl_graph.m @@ -0,0 +1,18 @@ +function gdl = bnet_to_gdl_graph(bnet) +% BNET_TO_GDL_GRAPH Convert a Bayesian network to a GDL graph +% gdl = bnet_to_gdl_graph(bnet) +% +% Each node in the BN gets converted to a single node in the GDL graph, +% representing its family; its kernel function is the corresponding CPD. + +N = length(bnet.dag); +doms = cell(1,N); +for i=1:N + doms{i} = family(bnet.dag, i); +end + +U = mk_undirected(bnet.dag); +gdl = mk_gdl_graph(U, doms, bnet.node_sizes, bnet.CPD, 'equiv_class', bnet.equiv_class, ... + 'discrete', bnet.dnodes, 'chance', bnet.chance_nodes, ... + 'decision', bnet.decision_nodes, 'utility', bnet.utility_nodes); + diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe.m new file mode 100644 index 00000000..5f55e708 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe.m @@ -0,0 +1,58 @@ +function [mpe, ll] = calc_mpe(engine, evidence, break_ties) +% CALC_MPE Computes the most probable explanation of the evidence +% [mpe, ll] = calc_mpe_given_inf_engine(engine, evidence, break_ties) +% +% INPUT +% engine must support max-propagation +% evidence{i} is the observed value of node i, or [] if hidden +% break_ties is optional. If 1, we will force ties to be broken consistently +% by calling enter_evidence N times. +% +% OUTPUT +% mpe{i} is the most likely value of node i (cell array!) +% ll is the log-likelihood of the globally best assignment +% +% This currently only works when all hidden nodes are discrete + +if nargin < 3, break_ties = 0; end + + +[engine, ll] = enter_evidence(engine, evidence, 'maximize', 1); + +observed = ~isemptycell(evidence); + +if 0 % fgraphs don't support bnet_from_engine +onodes = find(observed); +bnet = bnet_from_engine(engine); +pot_type = determine_pot_type(bnet, onodes); +assert(pot_type == 'd'); +end + +scalar = 1; +evidence = evidence(:); % hack to handle unrolled DBNs +N = length(evidence); +mpe = cell(1,N); +for i=1:N + m = marginal_nodes(engine, i); + % observed nodes are all set to 1 inside the inference engine, so we must undo this + if observed(i) + mpe{i} = evidence{i}; + else + mpe{i} = argmax(m.T); + % Bug fix by Ron Zohar, 8/15/01 + % If there are ties, we must break them as follows (see Jensen96, p106) + if break_ties + evidence{i} = mpe{i}; + [engine, ll] = enter_evidence(engine, evidence, 'maximize', 1); + end + end + if length(mpe{i}) > 1, scalar = 0; end +end + +if nargout >= 2 + bnet = bnet_from_engine(engine); + ll = log_lik_complete(bnet, mpe(:)); +end +if 0 % scalar + mpe = cell2num(mpe); +end diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_bucket.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_bucket.m new file mode 100644 index 00000000..40602725 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_bucket.m @@ -0,0 +1,160 @@ +function [mpe, ll] = calc_mpe_bucket(bnet, new_evidence, max_over) +% +% PURPOSE: +% CALC_MPE Computes the most probable explanation to the network nodes +% given the evidence. +% +% [mpe, ll] = calc_mpe(engine, new_evidence, max_over) +% +% INPUT: +% bnet - the bayesian network +% new_evidence - optional, if specified - evidence to be incorporated [cell(1,n)] +% max_over - optional, if specified determines the variable elimination order [1:n] +% +% OUTPUT: +% mpe - the MPE assignmet for the net variables (or [] if no satisfying assignment) +% ll - log assignment probability. +% +% Notes: +% 1. Adapted from '@var_elim_inf_engine\marginal_nodes' for MPE by Ron Zohar, 8/7/01 +% 2. Only discrete potentials are supported at this time. +% 3. Complexity: O(nw*) where n is the number of nodes and w* is the induced tree width. +% 4. Implementation based on: +% - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference", +% UA1 96, pp. 211-219. + + +ns = bnet.node_sizes; +n = length(bnet.dag); +evidence = cell(1,n); +if (nargin<2) + new_evidence = evidence; +end + +onodes = find(~isemptycell(new_evidence)); % observed nodes +hnodes = find(isemptycell(new_evidence)); % hidden nodes +pot_type = determine_pot_type(bnet, onodes); + +if pot_type ~= 'd' + error('only disrete potentials supported at this time') +end + +for i=1:n + fam = family(bnet.dag, i); + CPT{i} = convert_to_pot(bnet.CPD{bnet.equiv_class(i)}, pot_type, fam(:), evidence); +end + +% handle observed nodes: set impossible cases' probability to zero +% rather than prun matrix (this makes backtracking easier) + +for ii=onodes + lIdx = 1:ns(ii); + lIdx = setdiff(lIdx, new_evidence{ii}); + + sCPT=struct(CPT{ii}); % violate object privacy + + sargs = ''; + for jj=1:(length(sCPT.domain)-1) + sargs = [sargs, ':,']; + end + for jj=lIdx + eval(['sCPT.T(', sargs, num2str(jj), ')=0;']); + end + CPT{ii}=dpot(sCPT.domain, sCPT.sizes, sCPT.T); +end + +B = cell(1,n); +for b=1:n + B{b} = mk_initial_pot(pot_type, [], [], [], []); +end + +if (nargin<3) + max_over = (1:n); +end +order = max_over; % no attempt to optimize this + + +% Initialize the buckets with the CPDs assigned to them +for i=1:n + b = bucket_num(domain_pot(CPT{i}), order); + B{b} = multiply_pots(B{b}, CPT{i}); +end + +% Do backward phase +max_over = max_over(length(max_over):-1:1); % reverse +for i=max_over(1:end-1) + % max-ing over variable i which occurs in bucket j + j = bucket_num(i, order); + rest = mysetdiff(domain_pot(B{j}), i); + %temp = marginalize_pot_max(B{j}, rest); + temp = marginalize_pot(B{j}, rest, 1); + b = bucket_num(domain_pot(temp), order); + % fprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b); + sB=struct(B{b}); % violate object privacy + if ~isempty(sB.domain) + B{b} = multiply_pots(B{b}, temp); + else + B{b} = temp; + end +end +result = B{1}; +marginal = pot_to_marginal(result); +[prob, mpe] = max(marginal.T); + +% handle impossible cases +if ~(prob>0) + mpe = []; + ll = -inf; + %warning('evidence has zero probability') + return +end + +ll = log(prob); + +% Do forward phase +for ii=2:n + marginal = pot_to_marginal(B{ii}); + mpeidx = []; + for jj=order(1:length(mpe)) + assert(ismember(jj, marginal.domain)) %%% bug + temp = find_equiv_posns(jj, marginal.domain); + mpeidx = [mpeidx, temp] ; + if isempty(temp) + mpeidx = [mpeidx, Inf] ; + end + end + [mpeidxsorted sortedtompe] = sort(mpeidx) ; + + % maximize the matrix obtained from assigning values from previous buckets. + % this is done by building a string and using eval. + + kk=1; + sargs = '('; + for jj=1:length(marginal.domain) + if (jj~=1) + sargs = [sargs, ',']; + end + if (mpeidxsorted(kk)==jj) + sargs = [sargs, num2str(mpe(sortedtompe(kk)))]; + if (kk<length(mpe)) + kk = kk+1 ; + end + else + sargs = [sargs, ':']; + end + end + sargs = [sargs, ')'] ; + eval(['[val, loc] = max(marginal.T', sargs, ');']) + mpe = [mpe loc]; +end +[I,J] = sort(order); +mpe = mpe(J); + + + +%%%%%%%%% + +function b = bucket_num(domain, order) + +b = max(find_equiv_posns(domain, order)); + diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_dbn.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_dbn.m new file mode 100644 index 00000000..8889d49f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_dbn.m @@ -0,0 +1,41 @@ +function [mpe, ll] = calc_mpe_dbn(engine, evidence, break_ties) +% CALC_MPE Computes the most probable explanation of the evidence +% [mpe, ll] = calc_mpe_dbn(engine, evidence, break_ties) +% +% INPUT +% engine must support max-propagation +% evidence{i,t} is the observed value of node i in slice t, or [] if hidden +% +% OUTPUT +% mpe{i,t} is the most likely value of node i (cell array!) +% ll is the log-likelihood of the globally best assignment +% +% This currently only works when all hidden nodes are discrete + +if nargin < 3, break_ties = 0; end + +if break_ties + disp('warning: break ties is ignored') +end + +[engine, ll] = enter_evidence(engine, evidence, 'maximize', 1); + +observed = ~isemptycell(evidence); +[ss T] = size(evidence); +scalar = 1; +N = length(evidence); +mpe = cell(ss,T); +bnet = bnet_from_engine(engine); +ns = bnet.node_sizes; +for t=1:T + for i=1:ss + m = marginal_nodes(engine, i, t); + % observed nodes are all set to 1 inside the inference engine, so we must undo this + if observed(i,t) + mpe{i,t} = evidence{i,t}; + else + assert(length(m.T) == ns(i)); + mpe{i,t} = argmax(m.T); + end + end +end diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_given_inf_engine.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_given_inf_engine.m new file mode 100644 index 00000000..cd17a623 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_given_inf_engine.m @@ -0,0 +1,32 @@ +function [mpe, prob] = calc_mpe_given_inf_engine(engine, evidence) +% CALC_MPE_GIVEN_ENGINE Computes the most probable explanation of the evidence +% [mpe, prob] = calc_mpe_given_inf_engine(engine, evidence) +% +% INPUT +% engine must support max-propagation +% evidence{i} is the obsevred value of node i, or [] if hidden +% +% OUTPUT +% mpe(i) is the most likely value of node i +% prob is the likelihood of the globally best assignment +% +% This currently only works when all nodes are discrete + +[engine, ll] = enter_evidence(engine, evidence); + +observed = ~isemptycell(evidence); +N = length(evidence); +mpe = zeros(1,N); +for i=1:N + m = marginal_nodes(engine, i); + % discrete observed nodes are all set to 1 inside the inference engine, so we must undo this + if observed(i) + mpe(i) = evidence{i}; + else + mpe(i) = argmax(m.T); + end +end + +bnet = bnet_from_engine(engine); +ll = log_lik_complete(bnet, num2cell(mpe(:))); +prob = exp(ll); diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_global.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_global.m new file mode 100644 index 00000000..cf191966 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_global.m @@ -0,0 +1,28 @@ +function [mpe, ll] = calc_mpe_global(bnet, evidence) +% CALC_MPE_GLOBAL Compute the most probable explanation(s) from the global joint +% [mpe, ll] = calc_mpe_global(bnet, evidence) +% +% mpe(k,i) is the most probable value of node i in the k'th global mode +% ll is the log likelihood +% +% We assume all nodes are discrete + +engine = global_joint_inf_engine(bnet); +engine = enter_evidence(engine, evidence); +S1 = struct(engine); % violate object privacy +S2 = struct(S1.jpot); % joint potential +prob = max(S2.T(:)); +modes = find(S2.T(:) == prob); + +ens = bnet.node_sizes; +onodes = find(~isemptycell(evidence)); +ens(onodes) = 1; +mpe = ind2subv(ens, modes); +for k=1:length(modes) + for i=onodes(:)' + mpe(k,i) = evidence{i}; + end +end +ll = log(prob); + +mpe = num2cell(mpe); diff --git a/sourcecodes/bnt-master/BNT/general/Old/compute_interface_nodes.m b/sourcecodes/bnt-master/BNT/general/Old/compute_interface_nodes.m new file mode 100644 index 00000000..85f5c693 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/compute_interface_nodes.m @@ -0,0 +1,31 @@ +function [int, persist, transient] = compute_interface_nodes(intra, inter) +% COMPUTE_INTERFACE_NODES Find the nodes in a DBN that represent a sufficient statistic +% [int, persist, transient] = compute_interface_nodes(intra, inter) +% +% The interface nodes are all those that has an incoming temporal arc, +% or which have a child which has an incoming temporal arc, +% where a temporal arc means one coming from the previous slice. +% (The parents of nodes with incoming temporal arcs are needed +% because moralization will bring them into the clique.) +% +% The persisent nodes are all those that have one or more incoming temporal arc. +% The transient nodes are all the non-persistent. +% +% See U. Kjaerulff, "dHugin: A computational system for dynamic +% time-sliced Bayesian networks", Intl. J. Forecasting (11) 89-111, 1995 + +n = length(intra); +int = []; +persist = []; +for u=1:n + if any(inter(:,u)) + int = [int u]; + persist = [persist u]; + end + if any(inter(:, children(intra, u))) + int = [int u]; + end +end +int = unique(int); +persist = unique(persist); +transient = mysetdiff(1:n, persist); diff --git a/sourcecodes/bnt-master/BNT/general/Old/mk_gdl_graph.m b/sourcecodes/bnt-master/BNT/general/Old/mk_gdl_graph.m new file mode 100644 index 00000000..ec5350d7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/Old/mk_gdl_graph.m @@ -0,0 +1,86 @@ +function gdl = mk_gdl_graph(G, domains, node_sizes, kernels, varargin) +% MK_GDL_GRAPH Make a GDL (generalized distributed law) graph +% gdl = mk_gdl_graph(G, domains, node_sizes, kernels, ...) +% +% A GDL graph is like a moralized, but untriangulated, Bayes net: +% each "node" represents a domain with a corresponding kernel function. +% For details, see "The Generalized Distributive Law", Aji and McEliece, +% IEEE Trans. Info. Theory, 46(2): 325--343, 2000 +% +% G(i,j) = 1 if there is an (undirected) edge between domains i,j +% +% domains{i} is the domain of node i +% +% node_sizes(i) is the number of values node i can take on, +% or the length of node i if i is a continuous-valued vector. +% node_sizes(i) = 1 if i is a utility node. +% +% kernels is the list of kernel functions +% +% The list below gives optional arguments [default value in brackets]. +% +% equiv_class - equiv_class(i)=j means factor node i gets its params from factors{j} [1:F] +% discrete - the list of nodes which are discrete random variables [1:N] +% chance - the list of nodes which are random variables [1:N] +% decision - the list of nodes which are decision nodes [ [] ] +% utility - the list of nodes which are utility nodes [ [] ] + + +ns = node_sizes; +N = length(domains); +vars = []; +for i=1:N + vars = myunion(vars, domains{i}); +end +Nvars = length(vars); + +gdl.equiv_class = 1:length(kernels); +gdl.chance_nodes = 1:Nvars; +gdl.utility_nodes = []; +gdl.decision_nodes = []; +gdl.dnodes = 1:Nvars; + +if nargin >= 5 + args = varargin; + nargs = length(args); + for i=1:2:nargs + switch args{i}, + case 'equiv_class', bnet.equiv_class = args{i+1}; + case 'chance', bnet.chance_nodes = args{i+1}; + case 'utility', bnet.utility_nodes = args{i+1}; + case 'decision', bnet.decision_nodes = args{i+1}; + case 'discrete', bnet.dnodes = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + + +gdl.G = G; +gdl.vars = vars; +gdl.doms = domains; +gdl.node_sizes = node_sizes; +gdl.cnodes = mysetdiff(vars, gdl.dnodes); +gdl.kernels = kernels; +gdl.type = 'gdl'; + +% Compute a bit vector representation of the set of domains +% dom_bitv(i,j) = 1 iff variable j occurs in domain i +gdl.dom_bitv = zeros(N, length(vars)); +for i=1:N + gdl.dom_bitv(i, domains{i}) = 1; +end + +% compute the interesection of the domains on either side of each edge (separating set) +gdl.sepset = cell(N, N); +gdl.nbrs = cell(1,N); +for i=1:N + nbrs = neighbors(G, i); + gdl.nbrs{i} = nbrs; + for j = nbrs(:)' + gdl.sepset{i,j} = myintersect(domains{i}, domains{j}); + end +end + + diff --git a/sourcecodes/bnt-master/BNT/general/add_ev_to_dmarginal.m b/sourcecodes/bnt-master/BNT/general/add_ev_to_dmarginal.m new file mode 100644 index 00000000..b51f4983 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/add_ev_to_dmarginal.m @@ -0,0 +1,15 @@ +function fmarginal = add_ev_to_dmarginal(fmarginal, evidence, ns) +% ADD_EV_TO_DMARGINAL 'pump up' observed nodes back to their original size. +% fmarginal = add_ev_to_dmarginal(fmarginal, evidence, ns) +% +% We introduce 0s into the array in positions which are incompatible with the evidence. + +dom = fmarginal.domain; +odom = dom(~isemptycell(evidence(dom))); +vals = cat(1, evidence{odom}); +index = mk_multi_index(length(dom), find_equiv_posns(odom, dom), vals); +T = 0*myones(ns(dom)); +ens = ns(:)'; +ens(odom) = 1; +T(index{:}) = myreshape(fmarginal.T, ens(dom)); +fmarginal.T = T; diff --git a/sourcecodes/bnt-master/BNT/general/add_evidence_to_gmarginal.m b/sourcecodes/bnt-master/BNT/general/add_evidence_to_gmarginal.m new file mode 100644 index 00000000..5cbd7d96 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/add_evidence_to_gmarginal.m @@ -0,0 +1,78 @@ +function fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes) +% ADD_EVIDENCE_TO_GMARGINAL 'pump up' observed nodes back to their original size. +% function fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes) +% +% We introduce 0s into the array in positions which are incompatible with the evidence. +% for both discrete and continuous nodes. +% +% See also add_ev_to_dmarginal + +dom = fmarginal.domain; +fullm.domain = fmarginal.domain; + +% Find out which values of the discrete parents (if any) are compatible with +% the discrete evidence (if any). +dnodes = mysetdiff(1:length(ns), cnodes); +ddom = myintersect(dom, dnodes); +cdom = myintersect(dom, cnodes); +odom = dom(~isemptycell(evidence(dom))); +hdom = dom(isemptycell(evidence(dom))); + +% Find the entries in the big table that are compatible with the discrete evidence. +% (We will put the probabilities from the small inferred table into these positions.) +% We could use add_ev_to_dmarginal to do this. +dobs = myintersect(ddom, odom); +dvals = cat(1, evidence{dobs}); +ens = ns; % effective node sizes +ens(dobs) = 1; +S = prod(ens(ddom)); +subs = ind2subv(ens(ddom), 1:S); +mask = find_equiv_posns(dobs, ddom); +%subs(mask) = dvals; % bug fix by P. Brutti +for i=1:length(mask), + subs(:,mask(i)) = dvals(i); +end +supportedQs = subv2ind(ns(ddom), subs); + +if isempty(ddom) + Qarity = 1; +else + Qarity = prod(ns(ddom)); +end +fullm.T = zeros(Qarity, 1); +fullm.T(supportedQs) = fmarginal.T(:); +fullm.T = myreshape(fullm.T, ns(ddom)); + + +if isempty(cdom) + fullm.mu = []; + fullm.sigma = []; + return; +end + +% Now put the hidden cts parts into their right blocks, +% leaving the observed cts parts as 0. +cobs = myintersect(cdom, odom); +chid = myintersect(cdom, hdom); +cvals = cat(1, evidence{cobs}); +n = sum(ns(cdom)); +fullm.mu = zeros(n,Qarity); +fullm.Sigma = zeros(n,n,Qarity); + +if ~isempty(chid) + chid_blocks = block(find_equiv_posns(chid, cdom), ns(cdom)); +end +if ~isempty(cobs) + cobs_blocks = block(find_equiv_posns(cobs, cdom), ns(cdom)); +end + +for i=1:length(supportedQs) + Q = supportedQs(i); + if ~isempty(chid) + fullm.mu(chid_blocks, Q) = fmarginal.mu(:, i); + fullm.Sigma(chid_blocks, chid_blocks, Q) = fmarginal.Sigma(:,:,i); + end + if ~isempty(cobs) + fullm.mu(cobs_blocks, Q) = cvals(:); + end +end diff --git a/sourcecodes/bnt-master/BNT/general/bnet_to_fgraph.m b/sourcecodes/bnt-master/BNT/general/bnet_to_fgraph.m new file mode 100644 index 00000000..fa2763b9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/bnet_to_fgraph.m @@ -0,0 +1,16 @@ +function fg = bnet_to_fgraph(bnet) +% BNET_TO_FGRAPH Convert a Bayes net to a factor graph +% fg = bnet_to_fgraph(bnet) +% +% We create one factor per family, whose kernel is the CPD + +nnodes = length(bnet.dag); +G = zeros(nnodes, nnodes); +for i=1:nnodes + G(family(bnet.dag, i), i) = 1; +end + +fg = mk_fgraph(G, bnet.node_sizes, bnet.CPD, 'equiv_class', bnet.equiv_class, 'discrete', bnet.dnodes); + + + diff --git a/sourcecodes/bnt-master/BNT/general/compute_fwd_interface.m b/sourcecodes/bnt-master/BNT/general/compute_fwd_interface.m new file mode 100644 index 00000000..db4807c7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/compute_fwd_interface.m @@ -0,0 +1,11 @@ +function int = compute_fwd_interface(intra, inter) +% COMPUTE_FWD_INTERFACE Compute nodes with children in the next slice +% function int = compute_fwd_interface(intra, inter) + +int = []; +ss = length(intra); +for u=1:ss + if any(inter(u,:)) + int = [int u]; + end +end diff --git a/sourcecodes/bnt-master/BNT/general/compute_interface_nodes.m b/sourcecodes/bnt-master/BNT/general/compute_interface_nodes.m new file mode 100644 index 00000000..43370172 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/compute_interface_nodes.m @@ -0,0 +1,42 @@ +function [interface, persist, transient] = compute_interface_nodes(intra, inter) +% COMPUTE_INTERFACE_NODES Find the nodes in a DBN that represent a sufficient statistic +% [interface, persist, transient] = compute_interface_nodes(intra, inter) +% +% The interface nodes are all those that has an incoming temporal arc, +% or which are parents of such nodes. +% If the parents are in the previous slice, this just means they have an +% outgoing temporal arc. +% (The parents of nodes with incoming temporal arcs are needed +% because moralization will bring them into the clique.) +% +% The persisent nodes are all those that have one or more incoming temporal arc. +% The transient nodes are all the non-persistent. +% +% See U. Kjaerulff, "dHugin: A computational system for dynamic +% time-sliced Bayesian networks", Intl. J. Forecasting (11) 89-111, 1995 + +n = length(intra); +interface = []; +persist = []; +% any nodes with incoming arcs +for u=1:n + if any(inter(:,u)) + interface = [interface u]; + persist = [persist u]; + end +end +% Any nodes which are parents of nodes with incoming arcs +for u=1:n + cs = children(intra, u); + if any(inter(:, cs)) + interface = [interface u]; + end + %cs = children(inter, u); + % if ~isempty(myintersect(cs, persist)) + % interface = [interface u]; + %end +end +interface = unique(interface); +persist = unique(persist); +transient = mysetdiff(1:n, persist); + diff --git a/sourcecodes/bnt-master/BNT/general/compute_joint_pot.m b/sourcecodes/bnt-master/BNT/general/compute_joint_pot.m new file mode 100644 index 00000000..056aa998 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/compute_joint_pot.m @@ -0,0 +1,17 @@ +function [jpot, loglik] = compute_joint_pot(bnet, nodes, evidence, domain) +% COMPUTE_JOINT_POT Compute the global joint potential of a Bayes net +% function jpot = compute_joint_pot(bnet, nodes, evidence, domain) + +if nargin < 4, domain = nodes; end + +onodes = find(~isemptycell(evidence)); +pot_type = determine_pot_type(bnet, onodes, domain); + +jpot = mk_initial_pot(pot_type, domain, bnet.node_sizes, bnet.cnodes, onodes); +for i=nodes(:)' + e = bnet.equiv_class(i); + fam = family(bnet.dag, i); + pot = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence); + jpot = multiply_by_pot(jpot, pot); +end +%[jpot, loglik] = normalize_pot(jpot); % causes errors in asia_dt1 etc diff --git a/sourcecodes/bnt-master/BNT/general/compute_minimal_interface.m b/sourcecodes/bnt-master/BNT/general/compute_minimal_interface.m new file mode 100644 index 00000000..34dc8e62 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/compute_minimal_interface.m @@ -0,0 +1,25 @@ +function clqs = compute_minimal_interface(intra, inter) + +int = compute_fwd_interface(intra, inter); +ss = length(intra); +Z = zeros(ss); +dag = [intra inter; + Z intra]; +G = moralize(dag); +intra2 = G(1:ss,1:ss); +inter2 = G(1:ss,(1:ss)+ss); +G = unroll_dbn_topology(intra2, inter2, ss); +T = ss; +last_slice = (1:ss) + (T-1)*ss; +G = (G + G')/2; % mk symmetric +G2 = (expm(full(G)) > 0); % closure of graph +G3 = G2(last_slice, last_slice); +[c,v] = scc(G3); % connected components +ncomp = size(v,1); +clqs = cell(1,ncomp); +for i=1:ncomp + ndx = find(v(i,:)>0); + clqs{i} = v(i,ndx); + clqs{i} = myintersect(clqs{i}, int); +end + diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_pots.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_pots.m new file mode 100644 index 00000000..425e23b7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_pots.m @@ -0,0 +1,30 @@ +function CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type, softCPDpot) +% CONVERT_DBN_CPDS_TO_POTS Convert CPDs of (possibly instantiated) DBN nodes to potentials +% CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type, softCPDpot) +% +% CPDpot{n,t} is a potential containing P(n,t|pa(n,t), ev) +% softCPDpot{n,t} is a potential containing P(n,t|pa(n,t), ev) insted of using n's CPD + +[ss T] = size(evidence); + +if nargin < 4, softCPDpot = cell(ss,T); end +CPDpot = softCPDpot; + +% Convert CPDs of instantiated nodes to potential form +t = 1; +for n=1:ss + fam = family(bnet.dag, n); + e = bnet.equiv_class(n, 1); + if isempty(softCPDpot{n,t}) + CPDpot{n,t} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence(:,1)); + end +end +for n=1:ss + fam = family(bnet.dag, n, 2); + e = bnet.equiv_class(n, 2); + for t=2:T + if isempty(softCPDpot{n,t}) + CPDpot{n,t} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence(:,t-1:t)); + end + end +end diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables.m new file mode 100644 index 00000000..c10d58f7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables.m @@ -0,0 +1,201 @@ +function CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence) +% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables +% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence) +% +% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev) +% All hidden nodes are assumed to be discrete. +% We assume the observed nodes are the same in every slice. +% +% Evaluating the conditional likelihood of long evidence sequences can be very slow, +% so we take pains to vectorize where possible. + +[ss T] = size(evidence); +%obs_bitv = ~isemptycell(evidence(:)); +obs_bitv = zeros(1, 2*ss); +obs_bitv(bnet.observed) = 1; +obs_bitv(bnet.observed+ss) = 1; + +ns = bnet.node_sizes(:); +CPDpot = cell(ss,T); + +for n=1:ss + % slice 1 + t = 1; + ps = parents(bnet.dag, n); + e = bnet.equiv_class(n, 1); + if ~any(obs_bitv(ps)) + CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t}); + else + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1)); + end + +% special cases: c=child, p=parents, d=discrete, h=hidden, 1sl=1slice +% if c=h=1 then c=d=1, since hidden nodes must be discrete +% c=h c=d p=h p=d 1sl method +% --------------------------- +% 1 1 1 1 - replicate CPT +% - 1 - 1 - evaluate CPT on evidence * +% 0 1 1 1 1 dhmm +% 0 0 1 1 1 ghmm +% other loop +% +% * = any subset of the domain may be observed + +% Example where all of the special cases occur - a hierarchical HMM +% where the top layer (G) and leaves (Y) are observed and +% all nodes are discrete except Y. +% (O turns on if Y is an outlier) + +% G ---------> G +% | | +% v v +% S --------> S +% | | +% v v +% Y Y +% ^ ^ +% | | +% O O + +% Evaluating P(yt|St,Ot) is the ghmm case +% Evaluating P(St|S(t-1),gt) is the eval CPT case +% Evaluating P(gt|g(t-1) is the eval CPT case (hdom = []) +% Evaluating P(Ot) is the replicated CPT case + +% Cts parents (e.g., inputs) would require an additional special case for speed + + + % slices 2..T + [ss T] = size(evidence); + self = n+ss; + ps = parents(bnet.dag, self); + e = bnet.equiv_class(n, 2); + + if 1 + debug = 0; + hidden_child = ~obs_bitv(n); + discrete_child = myismember(n, bnet.dnodes); + hidden_ps = all(~obs_bitv(ps)); + discrete_ps = mysubset(ps, bnet.dnodes); + parents_in_same_slice = all(ps > ss); + + if hidden_child & discrete_child & hidden_ps & discrete_ps + CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug); + elseif discrete_child & discrete_ps + CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug); + elseif discrete_child & hidden_ps & discrete_ps & parents_in_same_slice + CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug); + elseif ~discrete_child & hidden_ps & discrete_ps & parents_in_same_slice + CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug); + else + if debug, fprintf('node %d, slow\n', n); end + for t=2:T + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t)); + end + end + end + + if 0 + for t=2:T + CPDpot2{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t)); + if ~approxeq(CPDpot{n,t}, CPDpot2{n,t}) + fprintf('CPDpot n=%d, t=%d\n',n,t); + keyboard + end + end + end + + +end + + + + +%%%%%%% +function CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug) + +[ss T] = size(evidence); +if debug, fprintf('node %d, repl\n', n); end +e = bnet.equiv_class(n, 2); +CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []); +CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]); + + + +%%%%%%% +function CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug) + +[ss T] = size(evidence); +self = n+ss; +ps = parents(bnet.dag, self); +e = bnet.equiv_class(n, 2); +ns = bnet.node_sizes(:); +% Example: given CPT(p1, p2, p3, p4, c), where p1,p3 are observed +% we create CPT([p2 p4 c], [p1 p3]). +% We then convert all observed p1,p3 into indices ndx +% and return CPT(:, ndx) +CPT = CPD_to_CPT(bnet.CPD{e}); +domain = [ps self]; +% if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2, +% odom = [3 8], hdom = 7 +odom_rel = find(obs_bitv(domain)); +hdom_rel = find(~obs_bitv(domain)); +odom = domain(odom_rel); +hdom = domain(hdom_rel); +if isempty(hdom) + CPT = CPT(:); +else + CPT = permute(CPT, [hdom_rel odom_rel]); + CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom))); +end +parents_in_same_slice = all(ps > ss); +if parents_in_same_slice + if debug, fprintf('node %d eval 1 slice\n', n); end + data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1 +else + if debug, fprintf('node %d eval 2 slice\n', n); end + % there's probably a way of vectorizing this... + data = zeros(length(odom), T-1); + for t=2:T + ev = evidence(:,t-1:t); + ev = ev(:); + ev2 = ev(odom); + data(:,t-1) = cat(1, ev2{:}); + %data(:,t-1) = cell2num(ev2); + end +end +ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t) +if isempty(hdom) + CPDpot(n,2:T) = num2cell(CPT(ndx)); % a cell array of floats +else + CPDpot(n,2:T) = num2cell(CPT(:, ndx), 1); % a cell array of column vectors +end + +%%%%%%% +function CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug) + +if debug, fprintf('node %d, dhmm\n', n); end +[ss T] = size(evidence); +self = n+ss; +ps = parents(bnet.dag, self); +e = bnet.equiv_class(n, 2); +ns = bnet.node_sizes(:); +CPT = CPD_to_CPT(bnet.CPD{e}); +CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents? +%obslik = mk_dhmm_obs_lik(cell2num(evidence(n,2:T)), CPT); +obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT); +CPDpot(n,2:T) = num2cell(obslik, 1); + + +%%%%%%% +function CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug) + +if debug, fprintf('node %d, ghmm\n', n); end +[ss T] = size(evidence); +e = bnet.equiv_class(n, 2); +S = struct(bnet.CPD{e}); +ev2 = cell2num(evidence(n,2:T)); +%obslik = mk_ghmm_obs_lik(ev2, S.mean, S.cov); +obslik = eval_pdf_cond_gauss(ev2, S.mean, S.cov); +CPDpot(n,2:T) = num2cell(obslik, 1); + diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables1.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables1.m new file mode 100644 index 00000000..c1b2c756 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables1.m @@ -0,0 +1,162 @@ +function CPDpot = convert_dbn_CPDs_to_tables1(bnet, evidence) +% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables +% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence) +% +% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev) +% All hidden nodes are assumed to be discrete +% We assume the observed nodes are the same in every slice +% +% Evaluating the conditional likelihood of the evidence can be very slow, +% so we take pains to vectorize where possible, i.e., we try to avoid +% calling convert_to_table + +[ss T] = size(evidence); +%obs_bitv = ~isemptycell(evidence(:)); +obs_bitv = zeros(1, 2*ss); +obs_bitv(bnet.observed) = 1; +obs_bitv(bnet.observed+ss) = 1; + +ns = bnet.node_sizes(:); +CPDpot = cell(ss,T); + +for n=1:ss + % slice 1 + t = 1; + ps = parents(bnet.dag, n); + e = bnet.equiv_class(n, 1); + if ~any(obs_bitv(ps)) + CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t}); + else + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1)); + end + + % slices 2..T + debug = 1; + if ~obs_bitv(n) + CPDpot = helper_hidden_child(bnet, evidence, n, CPDpot, obs_bitv, debug); + else + CPDpot = helper_obs_child(bnet, evidence, n, CPDpot, obs_bitv, debug); + end +end + +if 0 +CPDpot2 = convert_dbn_CPDs_to_tables_slow(bnet, evidence); +for t=1:T + for n=1:ss + if ~approxeq(CPDpot{n,t}, CPDpot2{n,t}) + fprintf('CPDpot n=%d, t=%d\n',n,t); + keyboard + end + end +end +end + + +% special cases: c=child, p=parents, d=discrete, h=hidden, 1=1slice +% if c=h=1 then c=d=1, since hidden nodes must be discrete +% c=h c=d p=h p=d p=1 method +% --------------------------- +% 1 1 1 1 - replicate CPT +% 0 1 1 1 1 dhmm +% 0 0 1 1 1 ghmm +% - 1 - 1 - evaluate CPT on evidence +% other loop + +%%%%%%% +function CPDpot = helper_hidden_child(bnet, evidence, n, CPDpot, obs_bitv, debug) + +[ss T] = size(evidence); +self = n+ss; +ps = parents(bnet.dag, self); +e = bnet.equiv_class(n, 2); +ns = bnet.node_sizes(:); +if ~any(obs_bitv(ps)) % all parents are hidden (hence discrete) + if debug, fprintf('node %d is hidden, all ps are hidden\n', n); end + if myismember(n, bnet.dnodes) + %CPT = CPD_to_CPT(bnet.CPD{e}); + %CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); + CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []); + CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]); + else + error(['hidden cts node disallowed']) + end +else % some parents are observed - slow + if mysubset(ps, bnet.dnodes) % all parents are discrete + % given CPT(p1, p2, p3, p4, c), where p1,p3 are observed + % we create CPT([p2 p4 c], [p1 p3]). + % We then convert all observed p1,p3 into indices ndx + % and return CPT(:, ndx) + CPT = CPD_to_CPT(bnet.CPD{e}); + domain = [ps self]; + % if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2, + % odom = [3 8], hdom = 7 + odom_rel = find(obs_bitv(domain)); + hdom_rel = find(~obs_bitv(domain)); + odom = domain(odom_rel); + hdom = domain(hdom_rel); + CPT = permute(CPT, [hdom_rel odom_rel]); + CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom))); + parents_in_same_slice = all(ps > ss); + if parents_in_same_slice + if debug, fprintf('node %d is hidden, some ps are obs, all ps discrete, 1 slice\n', n); end + data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1 + else + if debug, fprintf('node %d is hidden, some ps are obs, all ps discrete, 2 slice\n', n); end + data = zeros(length(odom), T-1); + for t=2:T + ev = evidence(:,t-1:t); + data(:,t-1) = cell2num(ev(odom)); + end + end + ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t) + CPDpot(n,2:T) = num2cell(CPT(:, ndx), [1 2]); + else % some parents are cts - v slow + if debug, fprintf('node %d is hidden, some ps are obs, some ps cts\n', n); end + for t=2:T + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t)); + end + end +end + +%%%%%%% +function CPDpot = helper_obs_child(bnet, evidence, n, CPDpot, obs_bitv, debug) + +[ss T] = size(evidence); +self = n+ss; +ps = parents(bnet.dag, self); +e = bnet.equiv_class(n, 2); +ns = bnet.node_sizes(:); +if ~any(obs_bitv(ps)) % all parents are hidden + parents_in_same_slice = all(ps > ss); + if parents_in_same_slice + if debug, fprintf('node %d is obs, all ps are hidden, 1 slice\n', n); end + ps1 = ps - ss; + if myismember(n, bnet.dnodes) + CPT = CPD_to_CPT(bnet.CPD{e}); + CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents? + obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT); + CPDpot(n,2:T) = num2cell(obslik, 1); + else + S = struct(bnet.CPD{e}); + obslik = eval_pdf_cond_gauss(cell2num(evidence(n,2:T)), S.mean, S.cov); + CPDpot(n,2:T) = num2cell(obslik, 1); + end + else % parents span 2 slices - slow + if debug, fprintf('node %d is obs, all ps are hidden , 2 slice\n', n); end + for t=2:T + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t)); + end + end +else + if isempty(ps) % observed root + if debug, fprintf('node %d is obs, no ps\n', n); end + CPT = CPD_to_CPT(bnet.CPD{e}); + data = cell2num(evidence(n,2:T)); + CPDpot(n,2:T) = CPT(data); + else % some parents are observed - slow + if debug, fprintf('node %d is obs, some ps are obs\n', n); end + for t=2:T + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t)); + end + end +end diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables_slow.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables_slow.m new file mode 100644 index 00000000..90e87046 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables_slow.m @@ -0,0 +1,41 @@ +function CPDpot = convert_dbn_CPDs_to_tables_slow(bnet, evidence) +% CONVERT_DBN_CPDS_TO_TABLES_SLOW Convert CPDs of (possibly instantiated) DBN nodes to tables +% CPDpot = convert_dbn_CPDs_to_tables_slow(bnet, evidence) +% +% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev) +% All hidden nodes are assumed to be discrete +% +% Non-vectorized method; this is less efficient for long sequences of observed Gaussian +% nodes, because of the (unnecessary) repeated matrix inversion. + +obs_bitv = ~isemptycell(evidence(:)); +[ss T] = size(evidence); +ns = bnet.node_sizes(:); + +CPDpot = cell(ss,T); + +t = 1; +for n=1:ss + %ps = engine.bnet_parents{n}; + ps = parents(bnet.dag, n); + e = bnet.equiv_class(n, 1); + if ~any(obs_bitv(ps)) + CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t}); + else + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1)); + end +end +for t=2:T + for n=1:ss + self = n+ss; + ps = parents(bnet.dag, self); + e = bnet.equiv_class(n, 2); + if ~any(obs_bitv(ps)) + CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t}); + else + CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t)); + end + end +end + + diff --git a/sourcecodes/bnt-master/BNT/general/dbn_to_bnet.m b/sourcecodes/bnt-master/BNT/general/dbn_to_bnet.m new file mode 100644 index 00000000..e15b9309 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/dbn_to_bnet.m @@ -0,0 +1,13 @@ +function bnet = dbn_to_bnet(dbn, T) +% DBN_TO_BNET Convert a DBN to a static network by unroll for T slices +% bnet = dbn_to_bnet(dbn, T) + +ss = length(dbn.intra); +eclass = [dbn.equiv_class(:,1) repmat(dbn.equiv_class(:,2), 1, T-1)]; +dnodes = unroll_set(dbn.dnodes_slice, ss, T); +ns = repmat(dbn.node_sizes_slice(:), 1, T); +dag = unroll_dbn_topology(dbn.intra, dbn.inter, T, dbn.intra1); +onodes = unroll_set(dbn.observed(:), ss, T); +bnet = mk_bnet(dag, ns(:), 'discrete', dnodes(:), 'equiv_class', eclass(:), 'observed', onodes(:)); +bnet.CPD = dbn.CPD; + diff --git a/sourcecodes/bnt-master/BNT/general/dbn_to_hmm.m b/sourcecodes/bnt-master/BNT/general/dbn_to_hmm.m new file mode 100644 index 00000000..8a783999 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/dbn_to_hmm.m @@ -0,0 +1,81 @@ +function [startprob, transprob, obsprob] = dbn_to_hmm(bnet) +% DBN_TO_HMM % Convert DBN params to HMM params +% [startprob, transprob, obsprob] = dbn_to_hmm(bnet, onodes) +% startprob(i) +% transprob(i,j) +% obsprob{k}.big_CPT(i,o) if k'th observed node is discrete +% obsprob{k}.big_mu(:,i), .big_Sigma(:,:,i) if k'th observed node is Gaussian +% Big means the domain contains all the hidden discrete nodes, not just the parents. + +% Called by constructor and by update_engine + +ss = length(bnet.intra); +onodes = bnet.observed; +hnodes = mysetdiff(1:ss, onodes); +evidence = cell(ss, 2); +ns = bnet.node_sizes(:); +Qh = prod(ns(hnodes)); +tmp = dpot_to_table(compute_joint_pot(bnet, hnodes, evidence)); +startprob = reshape(tmp, Qh, 1); + +tmp = dpot_to_table(compute_joint_pot(bnet, hnodes+ss, evidence, [hnodes hnodes+ss])); +transprob = mk_stochastic(reshape(tmp, Qh, Qh)); + +% P(o|ps) is used by mk_hmm_obs_lik_vec for a single time slice +% P(o|h) (the big version), where h = all hidden nodes, is used by enter_evidence + +obsprob = cell(1, length(onodes)); +for i=1:length(onodes) + o = onodes(i); + if bnet.auto_regressive(o) + % We assume the parents of this node are all the hidden nodes in the slice, + % so the params already are "big". Also, we assume we regress only on our old selves. + % slice 1 + e = bnet.equiv_class(o); + CPD = struct(bnet.CPD{e}); + O = ns(o); + ps = bnet.parents{o}; + Qps = prod(ns(ps)); + obsprob{i}.big_mu0 = reshape(CPD.mean, [O Qps]); + obsprob{i}.big_Sigma0 = reshape(CPD.cov, [O O Qps]); + + % slice t>1 + e = bnet.equiv_class(o+ss); + CPD = struct(bnet.CPD{e}); + O = ns(o); + dps = mysetdiff(bnet.parents{o+ss}, o); + Qdps = prod(ns(dps)); + obsprob{i}.big_mu = reshape(CPD.mean, [O Qdps]); + obsprob{i}.big_Sigma = reshape(CPD.cov, [O O Qdps]); + obsprob{i}.big_W = reshape(CPD.weights, [O O Qdps]); + else + e = bnet.equiv_class(o+ss); + CPD = struct(bnet.CPD{e}); + O = ns(o); + ps = bnet.parents{o}; + Qps = prod(ns(ps)); + % We make a big potential, replicating the params if necessary + % e.g., for a 2 chain coupled HMM, mu(:,Q1) becomes mu(:,Q1,Q2) + bigpot = pot_to_marginal(compute_joint_pot(bnet, onodes(i), evidence, [hnodes onodes(i)])); + + if myismember(o, bnet.dnodes) + obsprob{i}.CPT = reshape(CPD.CPT, [Qps O]); + obsprob{i}.big_CPT = reshape(bigpot.T, Qh, O); + else + obsprob{i}.big_mu = bigpot.mu; + obsprob{i}.big_Sigma = bigpot.Sigma; + + if 1 + obsprob{i}.mu = reshape(CPD.mean, [O Qps]); + C = reshape(CPD.cov, [O O Qps]); + obsprob{i}.Sigma = C; + d = size(obsprob{i}.mu, 1); + for j=1:Qps + obsprob{i}.inv_Sigma(:,:,j) = inv(C(:,:,j)); + obsprob{i}.denom(j) = (2*pi)^(d/2)*sqrt(abs(det(C(:,:,j)))); + end + end + + end % if discrete + end % if ar +end % for diff --git a/sourcecodes/bnt-master/BNT/general/determine_elim_constraints.m b/sourcecodes/bnt-master/BNT/general/determine_elim_constraints.m new file mode 100644 index 00000000..01c3f14c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/determine_elim_constraints.m @@ -0,0 +1,43 @@ +function partial_order = determine_elim_constraints(bnet, onodes) +% DETERMINE_ELIM_CONSTRAINTS Determine what the constraints are (if any) on the elimination ordering. +% partial_order = determine_elim_constraints(bnet, onodes) +% +% A graph with different kinds of nodes (e.g., discrete and cts, or decision and rnd) is called marked. +% A strong root is guaranteed to exist if the marked graph is triangulated and does not have any paths of +% the form discrete -> cts -> discrete. In general we need to add extra edges to +% the moral graph to ensure this (see example in Lauritzen (1992) fig 3b). +% However, a simpler sufficient condition is to eliminate all the cts nodes before the discrete ones, +% because then, as we move from the leaves to the root, the cts nodes get marginalized away +% and we are left with purely discrete cliques. +% +% partial_order(i,j)=1 if we must marginalize j *before* i +% (so i will be nearer the strong root). +% If the hidden nodes are either all discrete or all cts, we set partial_order = []. +% +% For details, see +% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94. +% - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical +% association models", JASA 87(420):1098--1108, 1992. +% - K. Olesen, "Causal probabilistic networks with both discrete and continuous variables", +% IEEE Pami 15(3), 1993 + + +n = length(bnet.dag); +pot_type = determine_pot_type(bnet, onodes); +if (pot_type == 'd') || (pot_type == 'g') + partial_order = []; + return; +end + + +partial_order = sparse(n,n); +partial_order(bnet.dnodes, bnet.cnodes) = 1; + +% Integrate out cts nodes before their discrete parents - see Olesen (1993) p9 +% This method gives the wrong results on cg1.m! +if 0 +for i=bnet.cnodes(:)' + dps = myintersect(parents(bnet.dag, i), bnet.dnodes); + partial_order(dps, i)=1; +end +end diff --git a/sourcecodes/bnt-master/BNT/general/dispcpt.m b/sourcecodes/bnt-master/BNT/general/dispcpt.m new file mode 100644 index 00000000..ba4c5607 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/dispcpt.m @@ -0,0 +1,16 @@ +function display_CPT(CPT) + +n = ndims(CPT); +parents_size = size(CPT); +parents_size = parents_size(1:end-1); +child_size = size(CPT,n); +c = 1; +for i=1:prod(parents_size) + parent_inst = ind2subv(parents_size, i); + fprintf(1, '%d ', parent_inst); + fprintf(1, ': '); + index = num2cell([parent_inst 1]); + index{n} = ':'; + fprintf(1, '%6.4f ', CPT(index{:})); + fprintf(1, '\n'); +end diff --git a/sourcecodes/bnt-master/BNT/general/do_intervention.m b/sourcecodes/bnt-master/BNT/general/do_intervention.m new file mode 100644 index 00000000..2cfab871 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/do_intervention.m @@ -0,0 +1,13 @@ +function bnet = mutilate_bnet(bnet, nodes, vals) +% MUTILATE_BNET Clamp nodes to specific values (perform a surgical intervention) +% bnet = mutilate_bnet(bnet, nodes, vals) +% +% We make all the clamped nodes roots. + +ns = bnet.node_sizes; +for i=1:length(nodes) + X = nodes(i); + x = vals(i); + bnet.dag(:,X) = 0; + bnet.CPD{X} = root_CPD(bnet, X, x); +end diff --git a/sourcecodes/bnt-master/BNT/general/dsep.m b/sourcecodes/bnt-master/BNT/general/dsep.m new file mode 100644 index 00000000..d4d9c441 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/dsep.m @@ -0,0 +1,15 @@ +function sep = dsep(X, Y, S, G) +% DSEP Is X indep Y given S wrt DAG G? +% sep = dsep(X, Y, S, G) +% +% Instead of using the Bayes-Ball criterion, we see if S separates X and Y +% in the moralized ancestral graph. + +conn = reachability_graph(G); +M = myunion(myunion(X, Y), S); +[A,junk] = find(conn(:, M)); +A = unique(A); +A = myunion(A, M); +GM = moralize(G(A,A)); +%sep = graph_separated(GM, X, Y, S); +sep = graph_separated(GM, find_equiv_posns(X,A), find_equiv_posns(Y,A), find_equiv_posns(S,A)); diff --git a/sourcecodes/bnt-master/BNT/general/dsep_test.m b/sourcecodes/bnt-master/BNT/general/dsep_test.m new file mode 100644 index 00000000..f0f14638 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/dsep_test.m @@ -0,0 +1,15 @@ + +% Cowell et al p72 +G = zeros(10); +G(1,2)=1; +G(2,3)=1; +G(3,7)=1; +G(4,[5 8])=1; +G(5,6)=1; +G(6,7)=1; +G(7,[9 10])=1; +G(8,9)=1; + +dsep(1, 4, [5 7], G) +dsep(1, 4, [7], G) +dsep(1, 4, [10 5], G) diff --git a/sourcecodes/bnt-master/BNT/general/enumerate_scenarios.m b/sourcecodes/bnt-master/BNT/general/enumerate_scenarios.m new file mode 100644 index 00000000..7d143065 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/enumerate_scenarios.m @@ -0,0 +1,21 @@ +function [scenarios, log_probs] = enumerate_scenarios(bnet, evidence) +% ENUMERATE_SCENARIOS Enumerate all assignments, and return the prob. of the non-zeros ones +% function [scenarios, log_probs] = enumerate_scenarios(bnet, evidence) + +assert(isempty(bnet.cnodes)); +n = length(bnet.dag); +observed = ~isemptycell(evidence); +vals = cat(1,evidence{observed}); +vals = vals(:)'; +ns = bnet.node_sizes; + +log_probs = []; +scenarios = []; +for i=1:prod(ns) + inst = ind2subv(ns, i); % i'th instantiation + if isempty(vals) | inst(observed) == vals % agrees with evidence + ll = log_lik_complete(bnet, num2cell(inst(:))); + log_probs = [log_probs ll]; + scenarios = [scenarios(:)' inst]; + end +end diff --git a/sourcecodes/bnt-master/BNT/general/fgraph_to_bnet.m b/sourcecodes/bnt-master/BNT/general/fgraph_to_bnet.m new file mode 100644 index 00000000..9aaf40ef --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/fgraph_to_bnet.m @@ -0,0 +1,30 @@ +function bnet = fgraph_to_bnet(fg) +% FGRAPH_TO_BNET Convert a factor graph to a Bayes net +% bnet = fgraph_to_bnet(fg) +% +% We assume all factors are tabular_CPD. +% We create 1 dummy observed node for every factor. + +N = fg.nvars + fg.nfactors; +vnodes = 1:fg.nvars; +fnodes = fg.nvars+1:N; +dag = zeros(N); +for x=1:fg.nvars + dag(x, fnodes(fg.dep{x})) = 1; +end +ns = [fg.node_sizes ones(1, fg.nfactors)]; +discrete = [fg.dnodes fnodes]; +bnet = mk_bnet(dag, ns, 'discrete', discrete); +for x=1:fg.nvars + bnet.CPD{x} = tabular_CPD(bnet, x, 'CPT', 'unif'); +end +ev = cell(1, fg.nvars); % no evidence +for i=1:fg.nfactors + f = fnodes(i); + e = fg.equiv_class(i); + pot = convert_to_pot(fg.factors{e}, 'd', fg.dom{i}, ev); + m = pot_to_marginal(pot); + bnet.CPD{f} = tabular_CPD(bnet, f, 'CPT', m.T); +end + + diff --git a/sourcecodes/bnt-master/BNT/general/hodbn_to_bnet.m b/sourcecodes/bnt-master/BNT/general/hodbn_to_bnet.m new file mode 100644 index 00000000..3026665d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/hodbn_to_bnet.m @@ -0,0 +1,21 @@ +function bnet = hodbn_to_bnet(dbn, T) +% DBN_TO_BNET Convert a DBN to a static network by unroll for T slices +% bnet = dbn_to_bnet(dbn, T) +ss = length(dbn.intra); +[row,order] = size(dbn.equiv_class); +eclass = []; +for i = 1:min(order,T) + eclass = [eclass ; dbn.equiv_class(:,i)]; +end +if T > order + eclass = [eclass ; repmat(dbn.equiv_class(:,order),T-order,1)]; +end + +dnodes = unroll_set(dbn.dnodes_slice, ss, T); +ns = repmat(dbn.node_sizes_slice(:), 1, T); +dag = unroll_higher_order_topology(dbn.intra, dbn.inter, T, dbn.intra1); +onodes = unroll_set(dbn.observed(:), ss, T); +bnet = mk_bnet(dag, ns(:), 'discrete', dnodes(:), 'equiv_class', eclass(:), 'observed', onodes(:)); +bnet.CPD = dbn.CPD; + + diff --git a/sourcecodes/bnt-master/BNT/general/is_mnet.m b/sourcecodes/bnt-master/BNT/general/is_mnet.m new file mode 100644 index 00000000..c6202e41 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/is_mnet.m @@ -0,0 +1,3 @@ +function m = is_mnet(model) + +m = isfield(model, 'markov_net'); diff --git a/sourcecodes/bnt-master/BNT/general/linear_gaussian_to_cpot.m b/sourcecodes/bnt-master/BNT/general/linear_gaussian_to_cpot.m new file mode 100644 index 00000000..f5a4251b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/linear_gaussian_to_cpot.m @@ -0,0 +1,43 @@ +function pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence) +% LINEAR_GAUSSIAN_TO_CPOT Convert a linear Gaussian CPD to a canonical potential. +% pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence) +% +% We include any cts evidence, but ignore any discrete evidence. +% (Use gaussian_CPD_params_given_dps to use discrete evidence to select mu, Sigma, W.) + +odom = domain(~isemptycell(evidence(domain))); +hdom = domain(isemptycell(evidence(domain))); +cobs = myintersect(cnodes, odom); +chid = myintersect(cnodes, hdom); +cvals = cat(1, evidence{cobs}); + +%[g,h,K] = gaussian_to_canonical(mu, Sigma, W); +Sinv = inv(Sigma); +g = -0.5*mu'*Sinv*mu + log(normal_coef(Sigma)); +if isempty(W) || (size(W,2)==0) % no cts parents + h = Sinv*mu; + K = Sinv; +else + h = [-W'*Sinv*mu; Sinv*mu]; + K = [W'*Sinv*W -W'*Sinv'; + -Sinv*W Sinv]; +end + +if ~isempty(cvals) + %[g, h, K] = enter_evidence_canonical(g, h, K, chid, cobs, cvals(:), ns); + [hx, hy, KXX, KXY, KYX, KYY] = partition_matrix_vec(h, K, chid, cobs, ns); + y = cvals(:); + g = g + hy'*y - 0.5*y'*KYY*y; + if length(hx)==0 % isempty(X) % i.e., we have instantiated everything away + h = []; + K = []; + else + h = hx - KXY*y; + K = KXX; + end +end + +ns(odom) = 0; +pot = cpot(domain, ns(domain), g, h, K); + + diff --git a/sourcecodes/bnt-master/BNT/general/log_lik_complete.m b/sourcecodes/bnt-master/BNT/general/log_lik_complete.m new file mode 100644 index 00000000..de0d9c51 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/log_lik_complete.m @@ -0,0 +1,30 @@ +function L = log_lik_complete(bnet, cases, clamped) +% LOG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m), theta_i) for a completely observed data set +% L = log_lik_complete(bnet, cases, clamped) +% +% If there is a missing data, you must use an inference engine. +% cases(i,m) is the value assigned to node i in case m. +% (If there are vector-valued nodes, cases should be a cell array.) +% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros) +% Clamped nodes contribute a factor of 1.0 to the likelihood. + +if iscell(cases), usecell = 1; else usecell = 0; end + +n = length(bnet.dag); +ncases = size(cases, 2); +if n ~= size(cases, 1) + error('data should be of size nnodes * ncases'); +end + +if nargin < 3, clamped = zeros(n,ncases); end + +L = 0; +for i=1:n + ps = parents(bnet.dag, i); + e = bnet.equiv_class(i); + u = find(clamped(i,:)==0); + ll = log_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u)); + if approxeq(exp(ll), 0), fprintf('node %d has very low likelihood\n'); end + L = L + ll; +end + diff --git a/sourcecodes/bnt-master/BNT/general/log_marg_lik_complete.m b/sourcecodes/bnt-master/BNT/general/log_marg_lik_complete.m new file mode 100644 index 00000000..8737b757 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/log_marg_lik_complete.m @@ -0,0 +1,40 @@ +function L = log_marg_lik_complete(bnet, cases, clamped) +% LOG_MARG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m)) for a completely observed data set +% L = log_marg_lik_complete(bnet, cases, clamped) +% +% This differs from log_lik_complete because we integrate out the parameters. +% If there is a missing data, you must use an inference engine. +% cases(i,m) is the value assigned to node i in case m. +% (If there are vector-valued nodes, cases should be a cell array.) +% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros) +% Clamped nodes contribute a factor of 1.0 to the likelihood. +% +% If there is a single case, clamped is a list of the clamped nodes, not a bit vector. + +if iscell(cases), usecell = 1; else usecell = 0; end + +n = length(bnet.dag); +ncases = size(cases, 2); +if n ~= size(cases, 1) + error('data should be of size nnodes * ncases'); +end + +if ncases == 1 + if nargin < 3, clamped = []; end + clamp_set = clamped; + clamped = zeros(n,1); + clamped(clamp_set) = 1; +else + if nargin < 3, clamped = zeros(n,ncases); end +end + +L = 0; +for i=1:n + ps = parents(bnet.dag, i); + e = bnet.equiv_class(i); + u = find(clamped(i,:)==0); + L = L + log_marg_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u)); +end + + + diff --git a/sourcecodes/bnt-master/BNT/general/mk_bnet.m b/sourcecodes/bnt-master/BNT/general/mk_bnet.m new file mode 100644 index 00000000..37560ac2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_bnet.m @@ -0,0 +1,93 @@ +function bnet = mk_bnet(dag, node_sizes, varargin) +% MK_BNET Make a Bayesian network. +% +% BNET = MK_BNET(DAG, NODE_SIZES, ...) makes a graphical model with an arc from i to j iff DAG(i,j) = 1. +% Thus DAG is the adjacency matrix for a directed acyclic graph. +% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary. +% +% node_sizes(i) is the number of values node i can take on, +% or the length of node i if i is a continuous-valued vector. +% node_sizes(i) = 1 if i is a utility node. +% +% Below are the names of optional arguments [and their default value in brackets]. +% Pass as 'PropertyName1', PropertyValue1, 'PropertyName2', PropertyValue2, ... +% +% discrete - the list of nodes which are discrete random variables [1:N] +% equiv_class - equiv_class(i)=j means node i gets its params from CPD{j} [1:N] +% observed - the list of nodes which will definitely be observed in every case [ [] ] +% 'names' - a cell array of strings to be associated with nodes 1:n [{}] +% This creates an associative array, so you write e.g. +% 'evidence(bnet.names{'bar'}) = 42' instead of 'evidence(2} = 42' +% assuming names = { 'foo', 'bar', ...}. +% +% e.g., bnet = mk_bnet(dag, ns, 'discrete', [1 3]) +% +% For backwards compatibility with BNT2, you can also specify the parameters in the following order +% bnet = mk_bnet(dag, node_sizes, discrete_nodes, equiv_class) + +n = length(dag); + +% default values for parameters +bnet.equiv_class = 1:n; +bnet.dnodes = 1:n; % discrete +bnet.observed = []; +bnet.names = {}; + +if nargin >= 3 + args = varargin; + nargs = length(args); + if ~ischar(args{1}) + if nargs >= 1, bnet.dnodes = args{1}; end + if nargs >= 2, bnet.equiv_class = args{2}; end + else + for i=1:2:nargs + switch args{i}, + case 'equiv_class', bnet.equiv_class = args{i+1}; + case 'discrete', bnet.dnodes = args{i+1}; + case 'observed', bnet.observed = args{i+1}; + case 'names', bnet.names = assocarray(args{i+1}, num2cell(1:n)); + otherwise, + error(['invalid argument name ' args{i}]); + end + end + end +end + +bnet.observed = sort(bnet.observed); % for comparing sets +bnet.hidden = mysetdiff(1:n, bnet.observed(:)'); +bnet.hidden_bitv = zeros(1,n); +bnet.hidden_bitv(bnet.hidden) = 1; +bnet.dag = dag; +bnet.node_sizes = node_sizes(:)'; + +bnet.cnodes = mysetdiff(1:n, bnet.dnodes); +% too many functions refer to cnodes to rename it to cts_nodes - +% We hope it won't be confused with chance nodes! + +bnet.parents = cell(1,n); +for i=1:n + bnet.parents{i} = parents(dag, i); +end + +E = max(bnet.equiv_class); +mem = cell(1,E); +for i=1:n + e = bnet.equiv_class(i); + mem{e} = [mem{e} i]; +end +bnet.members_of_equiv_class = mem; + +bnet.CPD = cell(1, E); + +bnet.rep_of_eclass = zeros(1,E); +for e=1:E + mems = bnet.members_of_equiv_class{e}; + bnet.rep_of_eclass(e) = mems(1); +end + +directed = 1; +if ~acyclic(dag,directed) + error('graph must be acyclic') +end + +bnet.order = topological_sort(bnet.dag); diff --git a/sourcecodes/bnt-master/BNT/general/mk_dbn.m b/sourcecodes/bnt-master/BNT/general/mk_dbn.m new file mode 100644 index 00000000..82859bc3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_dbn.m @@ -0,0 +1,133 @@ +function bnet = mk_dbn(intra, inter, node_sizes, varargin) +% MK_DBN Make a Dynamic Bayesian Network. +% +% BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs +% from i in slice t to j in slice t iff intra(i,j) = 1, and +% from i in slice t to j in slice t+1 iff inter(i,j) = 1, +% for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1. +% node_sizes(i) is the number of values node i can take on. +% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary. +% See also mk_bnet. +% +% Optional arguments [default in brackets] +% 'discrete' - list of discrete nodes [1:n] +% 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ] +% 'eclass1' - equiv class for slice 1 [1:n] +% 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1] +% equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j}, +% i.e., nodes i and j have tied parameters. +% 'intra1' - topology of first slice, if different from others +% 'names' - a cell array of strings to be associated with nodes 1:n [{}] +% This creates an associative array, so you write e.g. +% 'evidence(bnet.names{'bar'}) = 42' instead of 'evidence(2} = 42' +% assuming names = { 'foo', 'bar', ...}. +% +% For backwards compatibility with BNT2, arguments can also be specified as follows +% bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1) +% +% After calling this function, you must specify the parameters (conditional probability +% distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc. + + +n = length(intra); +ss = n; +bnet.nnodes_per_slice = ss; +bnet.intra = intra; +bnet.inter = inter; +bnet.intra1 = intra; +dag = zeros(2*n); +dag(1:n,1:n) = bnet.intra1; +dag(1:n,(1:n)+n) = bnet.inter; +dag((1:n)+n,(1:n)+n) = bnet.intra; +bnet.dag = dag; +bnet.names = {}; + +directed = 1; +if ~acyclic(dag,directed) + error('graph must be acyclic') +end + + +bnet.eclass1 = 1:n; +%bnet.eclass2 = (1:n)+n; +bnet.eclass2 = bnet.eclass1; +for i=1:ss + if isequal(parents(dag, i+ss), parents(dag, i)+ss) + %fprintf('%d has isomorphic parents, eclass %d\n', i, bnet.eclass2(i)) + else + bnet.eclass2(i) = max(bnet.eclass2) + 1; + %fprintf('%d has non isomorphic parents, eclass %d\n', i, bnet.eclass2(i)) + end +end + +dnodes = 1:n; +bnet.observed = []; + +if nargin >= 4 + args = varargin; + nargs = length(args); + if ~isstr(args{1}) + if nargs >= 1, dnodes = args{1}; end + if nargs >= 2, bnet.eclass1 = args{2}; end + if nargs >= 3, bnet.eclass2 = args{3}; end + if nargs >= 4, bnet.intra1 = args{4}; end + else + for i=1:2:nargs + switch args{i}, + case 'discrete', dnodes = args{i+1}; + case 'observed', bnet.observed = args{i+1}; + case 'eclass1', bnet.eclass1 = args{i+1}; + case 'eclass2', bnet.eclass2 = args{i+1}; + case 'intra1', bnet.intra1 = args{i+1}; + %case 'ar_hmm', bnet.ar_hmm = args{i+1}; % should check topology + case 'names', bnet.names = assocarray(args{i+1}, num2cell(1:n)); + otherwise, + error(['invalid argument name ' args{i}]); + end + end + end +end + + +bnet.observed = sort(bnet.observed); % for comparing sets +ns = node_sizes; +bnet.node_sizes_slice = ns(:)'; +bnet.node_sizes = [ns(:) ns(:)]; + +cnodes = mysetdiff(1:n, dnodes); +bnet.dnodes_slice = dnodes; +bnet.cnodes_slice = cnodes; +bnet.dnodes = [dnodes dnodes+n]; +bnet.cnodes = [cnodes cnodes+n]; + +bnet.equiv_class = [bnet.eclass1(:) bnet.eclass2(:)]; +bnet.CPD = cell(1,max(bnet.equiv_class(:))); +eclass = bnet.equiv_class(:); +E = max(eclass); +bnet.rep_of_eclass = zeros(1,E); +for e=1:E + mems = find(eclass==e); + bnet.rep_of_eclass(e) = mems(1); +end + +ss = n; +onodes = bnet.observed; +hnodes = mysetdiff(1:ss, onodes); +bnet.hidden_bitv = zeros(1,2*ss); +bnet.hidden_bitv(hnodes) = 1; +bnet.hidden_bitv(hnodes+ss) = 1; + +bnet.parents = cell(1, 2*ss); +for i=1:ss + bnet.parents{i} = parents(bnet.dag, i); + bnet.parents{i+ss} = parents(bnet.dag, i+ss); +end + +bnet.auto_regressive = zeros(1,ss); +% ar(i)=1 means (observed) node i depends on i in the previous slice +for o=bnet.observed(:)' + if any(bnet.parents{o+ss} <= ss) + bnet.auto_regressive(o) = 1; + end +end + diff --git a/sourcecodes/bnt-master/BNT/general/mk_fgraph.m b/sourcecodes/bnt-master/BNT/general/mk_fgraph.m new file mode 100644 index 00000000..8ac50338 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_fgraph.m @@ -0,0 +1,60 @@ +function fg = mk_fgraph(G, node_sizes, factors, varargin) +% MK_FGRAPH Make a factor graph +% fg = mk_fgraph(G, node_sizes, factors, ...) +% +% A factor graph is a bipartite graph, with one side containing variables, +% and the other containing functions of (subsets of) these variables. +% For details, see "Factor Graphs and the Sum-Product Algorithm", +% F. Kschischang and B. Frey and H-A. Loeliger, +% IEEE Trans. Info. Theory, 2001 +% +% G(i,j) = 1 if there is an arc from variable i to factor j +% +% node_sizes(i) is the number of values node i can take on, +% or the length of node i if i is a continuous-valued vector. +% +% 'factors' is the list of factors (kernel functions) +% +% The list below gives optional arguments [default value in brackets]. +% +% equiv_class - equiv_class(i)=j means factor node i gets its params from factors{j} [1:F] +% discrete - the list of nodes which are discrete random variables [1:N] +% +% e.g., fg = mk_fgraph(G, [2 2], {bnet.CPD{1},bnet.CPD{2}}, 'discrete', [1 2]) + +fg.G = G; +fg.node_sizes = node_sizes; +fg.factors = factors; +[fg.nvars fg.nfactors] = size(G); + +% default values for parameters +fg.equiv_class = 1:fg.nfactors; +fg.dnodes = 1:fg.nvars; + +if nargin >= 4 + args = varargin; + nargs = length(args); + for i=1:2:nargs + switch args{i}, + case 'equiv_class', fg.equiv_class = args{i+1}; + case 'discrete', fg.dnodes = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + +% so that determine_pot_type will work... +fg.utility_nodes = []; +%fg.decision_nodes = []; +%fg.chance_nodes = fg.nvars; + +fg.dom = cell(1, fg.nfactors); +for f=1:fg.nfactors + fg.dom{f} = find(G(:,f)); +end +fg.dep = cell(1, fg.nvars); +for x=1:fg.nvars + fg.dep{x} = find(G(x,:)); +end +fg.cnodes = mysetdiff(1:fg.nvars, fg.dnodes); diff --git a/sourcecodes/bnt-master/BNT/general/mk_fgraph_given_ev.m b/sourcecodes/bnt-master/BNT/general/mk_fgraph_given_ev.m new file mode 100644 index 00000000..5e65ed10 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_fgraph_given_ev.m @@ -0,0 +1,48 @@ +function fg = mk_fgraph_given_ev(G, node_sizes, factors, ev_CPD, evidence, varargin) +% MK_FGRAPH_GIVEN_EV Make a factor graph where each node has its own private evidence term +% fg = mk_fgraph(G, node_sizes, factors, ev_CPD, evidence, ...) +% +% G, node_sizes and factors are as in mk_fgraph, but they refer to the hidden nodes. +% ev_CPD{i} is a CPD for the i'th hidden node; this will be converted into a factor +% for node i using evidence{i}. +% We currently assume all hidden nodes are discrete, for simplicity. +% +% The list below gives optional arguments [default value in brackets]. +% +% equiv_class - equiv_class(i)=j means factor node i gets its params from factors{j} [1:F] +% ev_equiv_class - ev_equiv_class(i)=j means evidence node i gets its params from ev_CPD{j} [1:N] + + +N = length(node_sizes); +nfactors = length(factors); + +% default values for parameters +eclass = 1:nfactors; +ev_eclass = 1:N; + +if nargin >= 6 + args = varargin; + nargs = length(args); + for i=1:2:nargs + switch args{i}, + case 'equiv_class', eclass = args{i+1}; + case 'ev_equiv_class', ev_eclass = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + +pot_type = 'd'; +for x=1:N + ev = cell(1,2); % cell 1 is the hidden parent, cell 2 is the observed child + ev(2) = evidence(x); + dom = 1:2; + F = convert_to_pot(ev_CPD{ev_eclass(x)}, pot_type, dom(:), ev); + M = pot_to_marginal(F); + %factors{end+1} = tabular_CPD('self', 1, 'ps', [], 'sz', node_sizes(x), 'CPT', M.T); + factors{end+1} = mk_isolated_tabular_CPD(node_sizes(x), {'CPT', M.T}); +end + +E = max(eclass); +fg = mk_fgraph([G eye(N)], node_sizes, factors, 'equiv_class', [eclass E+1:E+N]); diff --git a/sourcecodes/bnt-master/BNT/general/mk_higher_order_dbn.m b/sourcecodes/bnt-master/BNT/general/mk_higher_order_dbn.m new file mode 100644 index 00000000..eb01636e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_higher_order_dbn.m @@ -0,0 +1,181 @@ +function bnet = mk_higher_order_dbn(intra, inter, node_sizes, varargin) +% MK_DBN Make a Dynamic Bayesian Network. +% +% BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs +% from i in slice t to j in slice t iff intra(i,j) = 1, and +% from i in slice t to j in slice t+1 iff inter(i,j) = 1, +% for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1. +% node_sizes(i) is the number of values node i can take on. +% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary. +% See also mk_bnet. +% +% Optional arguments [default in brackets] +% 'discrete' - list of discrete nodes [1:n] +% 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ] +% 'eclass1' - equiv class for slice 1 [1:n] +% 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1] +% equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j}, +% i.e., nodes i and j have tied parameters. +% 'intra1' - topology of first slice, if different from others +% 'names' - a cell array of strings to be associated with nodes 1:n [{}] +% This creates an associative array, so you write e.g. +% 'evidence(bnet.names{'bar'}) = 42' instead of 'evidence(2} = 42' +% assuming names = { 'foo', 'bar', ...}. +% +% For backwards compatibility with BNT2, arguments can also be specified as follows +% bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1) +% +% After calling this function, you must specify the parameters (conditional probability +% distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc. + + +n = length(intra); +ss = n; +bnet.nnodes_per_slice = ss; +bnet.intra = intra; +bnet.inter = inter; +bnet.intra1 = intra; + +% As this method is used to generate a higher order Markov Model +% also connect from time slice t - i -> t with i > 1 has to be +% taken into account. + +%inter should be a three dimensional array where inter(:,:,i) +%describes the connections from time-slice t - i to t. +[rows,columns,order] = size(inter); +assert(rows == n); +assert(columns == n); +dag = zeros((order + 1)*n); + +i = 0; +while i <= order + j = i; + while j <= order + if j == i + dag(1 + i*n:(i+1)*n,1+i*n:(i+1)*n) = intra; + else + dag(1+i*n:(i+1)*n,1+j*n:(j+1)*n) = inter(:,:,j - i); + end + j = j + 1; + end; + i = i + 1; +end; + +bnet.dag = dag; +bnet.names = {}; + +directed = 1; +if ~acyclic(dag,directed) + error('graph must be acyclic') +end + +% Calculation of the equivalence classes +bnet.eclass1 = 1:n; +bnet.eclass = zeros(order + 1,ss); +bnet.eclass(1,:) = 1:n; +for i = 1:order + bnet.eclass(i+1,:) = bnet.eclass(i,:); + for j = 1:ss + if(isequal(parents(dag,(i-1)*n+j)+ss,parents(dag,(i*n + j)))) + %fprintf('%d has isomorphic parents, eclass %d \n',j,bnet.eclass(i,j)) + else + bnet.eclass(i + 1,j) = max(bnet.eclass(i+1,:))+1; + %fprintf('%d has non isomorphic parents, eclass %d \n',j,bnet.eclass(i,j)) + end; + end; +end; +bnet.eclass1 = 1:n; + +% To be compatible with whe rest of the code +bnet.eclass2 = bnet.eclass(2,:); + +dnodes = 1:n; +bnet.observed = []; + +if nargin >= 4 + args = varargin; + nargs = length(args); + if ~isstr(args{1}) + if nargs >= 1 dnodes = args{1}; end + if nargs >= 2 bnet.eclass1 = args{2}; bnet.eclass(1,:) = args{2}; end + if nargs >= 3 bnet.eclass2 = args{3}; bnet.eclass(2,:) = args{2}; end + if nargs >= 4 bnet.intra1 = args{4}; end + else + for i=1:2:nargs + switch args{i}, + case 'discrete', dnodes = args{i+1}; + case 'observed', bnet.observed = args{i+1}; + case 'eclass1', bnet.eclass1 = args{i+1}; bnet.eclass(1,:) = args{i+1}; + case 'eclass2', bnet.eclass2 = args{i+1}; bnet.eclass(2,:) = args{i+1}; + case 'eclass', bnet.eclass = args{i+1}; + case 'intra1', bnet.intra1 = args{i+1}; + %case 'ar_hmm', bnet.ar_hmm = args{i+1}; % should check topology + case 'names', bnet.names = assocarray(args{i+1}, num2cell(1:n)); + otherwise, + error(['invalid argument name ' args{i}]); + end + end + end +end + +bnet.observed = sort(bnet.observed); % for comparing sets +ns = node_sizes; +bnet.node_sizes_slice = ns(:)'; +bnet.node_sizes = repmat(ns(:),1,order + 1); + +cnodes = mysetdiff(1:n, dnodes); +bnet.dnodes_slice = dnodes; +bnet.cnodes_slice = cnodes; +bnet.dnodes = dnodes; +bnet.cnodes = cnodes; +% To adapt the function to higher order Markov models include dnodes for more +% time slices +for i = 1:order + bnet.dnodes = [bnet.dnodes dnodes+i*n]; + bnet.cnodes = [bnet.cnodes cnodes+i*n]; +end + +% Generieren einer Matrix, deren i-te Spalte die Aequivalenzklassen +% der i-ten Zeitscheibe enthaelt. +bnet.equiv_class = [bnet.eclass(1,:)]'; +for i = 2:(order + 1) + bnet.equiv_class = [bnet.equiv_class bnet.eclass(i,:)']; +end + +bnet.CPD = cell(1,max(bnet.equiv_class(:))); + +ss = n; +onodes = bnet.observed; +hnodes = mysetdiff(1:ss, onodes); +bnet.hidden_bitv = zeros(1,(order + 1)*ss); +for i = 0:order + bnet.hidden_bitv(hnodes +i*ss) = 1; +end; + +bnet.parents = cell(1, (order + 1)*ss); +for i=1:(order + 1)*ss + bnet.parents{i} = parents(bnet.dag, i); +end + +bnet.auto_regressive = zeros(1,ss); +% ar(i)=1 means (observed) node i depends on i in the previous slice +for o=bnet.observed(:)' + if any(bnet.parents{o+ss} <= ss) + bnet.auto_regressive(o) = 1; + end +end + + + + + + + + + + + + + + + diff --git a/sourcecodes/bnt-master/BNT/general/mk_limid.m b/sourcecodes/bnt-master/BNT/general/mk_limid.m new file mode 100644 index 00000000..3c10ed81 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_limid.m @@ -0,0 +1,93 @@ +function bnet = mk_limid(dag, node_sizes, varargin) +% MK_LIMID Make a limited information influence diagram +% +% BNET = MK_LIMID(DAG, NODE_SIZES, ...) +% DAG is the adjacency matrix for a directed acyclic graph. +% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary. +% For decision nodes, the parents must explicitely include all nodes +% on which it can depends, in contrast to the implicit no-forgetting assumption of influence diagrams. +% (For details, see "Representing and solving decision problems with limited information", +% Lauritzen and Nilsson, Management Science, 2001.) +% +% node_sizes(i) is the number of values node i can take on, +% or the length of node i if i is a continuous-valued vector. +% node_sizes(i) = 1 if i is a utility node. +% +% The list below gives optional arguments [default value in brackets]. +% +% chance - the list of nodes which are random variables [1:N] +% decision - the list of nodes which are decision nodes [ [] ] +% utility - the list of nodes which are utility nodes [ [] ] +% equiv_class - equiv_class(i)=j means node i gets its params from CPD{j} [1:N] +% +% e.g., limid = mk_limid(dag, ns, 'chance', [1 3], 'utility', [2]) + +n = length(dag); + +% default values for parameters +bnet.chance_nodes = 1:n; +bnet.equiv_class = 1:n; +bnet.utility_nodes = []; +bnet.decision_nodes = []; +bnet.dnodes = 1:n; % discrete + +if nargin >= 3 + args = varargin; + nargs = length(args); + if ~isstr(args{1}) + if nargs >= 1, bnet.dnodes = args{1}; end + if nargs >= 2, bnet.equiv_class = args{2}; end + else + for i=1:2:nargs + switch args{i}, + case 'equiv_class', bnet.equiv_class = args{i+1}; + case 'chance', bnet.chance_nodes = args{i+1}; + case 'utility', bnet.utility_nodes = args{i+1}; + case 'decision', bnet.decision_nodes = args{i+1}; + case 'discrete', bnet.dnodes = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end + end +end + +bnet.limid = 1; + +bnet.dag = dag; +bnet.node_sizes = node_sizes(:)'; + +bnet.cnodes = mysetdiff(1:n, bnet.dnodes); +% too many functions refer to cnodes to rename it to cts_nodes - +% We hope it won't be confused with chance nodes! + +bnet.parents = cell(1,n); +for i=1:n + bnet.parents{i} = parents(dag, i); +end + +E = max(bnet.equiv_class); +mem = cell(1,E); +for i=1:n + e = bnet.equiv_class(i); + mem{e} = [mem{e} i]; +end +bnet.members_of_equiv_class = mem; + +bnet.CPD = cell(1, E); + +% for e=1:E +% i = bnet.members_of_equiv_class{e}(1); % pick arbitrary member +% switch type{e} +% case 'tabular', bnet.CPD{e} = tabular_CPD(bnet, i); +% case 'gaussian', bnet.CPD{e} = gaussian_CPD(bnet, i); +% otherwise, error(['unrecognized CPD type ' type{e}]); +% end +% end + +directed = 1; +if ~acyclic(dag,directed) + error('graph must be acyclic') +end + +bnet.order = topological_sort(bnet.dag); diff --git a/sourcecodes/bnt-master/BNT/general/mk_mnet.m b/sourcecodes/bnt-master/BNT/general/mk_mnet.m new file mode 100644 index 00000000..38d4cab3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_mnet.m @@ -0,0 +1,14 @@ +function mnet = mk_mnet(graph, node_sizes, cliques, potentials) +% MK_MNET Make a Markov network (Markov Random Field) +% +% mnet = mk_mnet(adj_mat, node_sizes, cliques, potentials) +% +% cliques{i} is a list of the nodes in clq i +% potentials{i} is a dpot object corresponding to the potential for clique i +% + +mnet.markov_net = 1; +mnet.graph = graph; +mnet.node_sizes = node_sizes; +mnet.user_cliques = cliques; +mnet.user_potentials = potentials; diff --git a/sourcecodes/bnt-master/BNT/general/mk_mrf2.m b/sourcecodes/bnt-master/BNT/general/mk_mrf2.m new file mode 100644 index 00000000..4d948864 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_mrf2.m @@ -0,0 +1,8 @@ +function mrf2 = mk_mrf2(adj_mat, pot) +% MK_MRF2 Make a Markov random field with pairwise potentials +% function mrf2 = mk_mrf2(adj_mat, pot) +% +% pot{i,j}(k1,k2) + +mrf2.adj_mat = adj_mat; +mrf2.pot = pot; diff --git a/sourcecodes/bnt-master/BNT/general/mk_mutilated_samples.m b/sourcecodes/bnt-master/BNT/general/mk_mutilated_samples.m new file mode 100644 index 00000000..21a91063 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_mutilated_samples.m @@ -0,0 +1,38 @@ +function [data, clamped] = mk_mutilated_samples(bnet, ncases, max_clamp, usecell) +% GEN_MUTILATED_SAMPLES Do random interventions and then draw random samples +% [data, clamped] = gen_mutilated_samples(bnet, ncases, max_clamp, usecell) +% +% At each step, we pick a random subset of size 0 .. max_clamp, and +% clamp these nodes to random values. +% +% data(i,m) is the value of node i in case m. +% clamped(i,m) = 1 if node i in case m was set by intervention. + +if nargin < 4, usecell = 1; end + +ns = bnet.node_sizes; +n = length(bnet.dag); +if usecell + data = cell(n, ncases); +else + data = zeros(n, ncases); +end +clamped = zeros(n, ncases); + +csubsets = subsets(1:n, max_clamp, 0); % includes the empty set +distrib_cset = normalise(ones(1, length(csubsets))); + +for m=1:ncases + cset = csubsets{sample_discrete(distrib_cset)}; + nvals = prod(ns(cset)); + distrib_cvals = normalise(ones(1, nvals)); + cvals = ind2subv(ns(cset), sample_discrete(distrib_cvals)); + mutilated_bnet = do_intervention(bnet, cset, cvals); + ev = sample_bnet(mutilated_bnet); + if usecell + data(:,m) = ev; + else + data(:,m) = cell2num(ev); + end + clamped(cset,m) = 1; +end diff --git a/sourcecodes/bnt-master/BNT/general/mk_named_CPT.m b/sourcecodes/bnt-master/BNT/general/mk_named_CPT.m new file mode 100644 index 00000000..5f39a981 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_named_CPT.m @@ -0,0 +1,52 @@ +function CPT2 = mk_named_CPT(family_names, names, dag, CPT1) +% MK_NAMED_CPT Permute the dimensions of a CPT so they agree with the internal numbering convention +% CPT2 = mk_named_CPT(family_names, names, dag, CPT1) +% +% This is best explained by example. +% Consider the following directed acyclic graph +% +% C +% / \ +% R S +% \ / +% W +% +% where all arcs point down. +% When we create the CPT for node W, we consider S as its first parent, and R as its +% second, and hence write +% +% S R W +% CPT1(1,1,:) = [1.0 0.0]; +% CPT1(2,1,:) = [0.2 0.8]; % P(W=1 | R=1, S=2) = 0.2 +% CPT1(1,2,:) = [0.1 0.9]; +% CPT1(2,2,:) = [0.01 0.99]; +% +% However, when we create the dag using mk_adj_mat, the nodes get topologically sorted, +% and by chance, node R preceeds node S in this ordering. +% Hence we should have written +% +% R S W +% CPT2(1,1,:) = [1.0 0.0]; +% CPT2(2,1,:) = [0.1 0.9]; +% CPT2(1,2,:) = [0.2 0.8]; % P(W=1 | R=1, S=2) = 0.2 +% CPT2(2,2,:) = [0.01 0.99]; +% +% Since we do not know the order of the nodes in advance, we can write +% CPT2 = mk_named_CPT({'S', 'R', 'W'}, names, dag, CPT1) +% where 'S', 'R', 'W' are the order of the dimensions we assumed (the child node must be last in this list), +% and names{i} is the name of the i'th node. + +n = length(family_names); +family_nums = zeros(1,n); +for i=1:n + family_nums(i) = stringmatch(family_names{i}, names); % was strmatch +end + +fam = family(dag, family_nums(end)); +perm = zeros(1,n); +for i=1:n + % perm(i) = find(family_nums(i) == fam); + perm(i) = find(fam(i) == family_nums); +end + +CPT2 = permute(CPT1, perm); diff --git a/sourcecodes/bnt-master/BNT/general/mk_slice_and_half_dbn.m b/sourcecodes/bnt-master/BNT/general/mk_slice_and_half_dbn.m new file mode 100644 index 00000000..80faad17 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/mk_slice_and_half_dbn.m @@ -0,0 +1,28 @@ +function bnet15 = mk_slice_and_half_dbn(bnet, int) +% function bnet = mk_slice_and_half_dbn(bnet, int) +% function bnet = mk_slice_and_half_dbn(bnet, int) +% +% Create a "1.5 slice" jtree, containing the interface nodes of slice 1 +% and all the nodes of slice 2 +% To keep the node numbering the same, we simply disconnect the non-interface nodes +% from slice 1, and set their size to 1. +% We do this to speed things up, and so that the likelihood is computed correctly. +% We do not need to do +% this if we just want to compute marginals (i.e., we can include nodes whose potentials will +% be left as all 1s). + +intra15 = bnet.intra; +ss = length(bnet.intra); +nonint = mysetdiff(1:ss, int); +for i=nonint(:)' + intra15(:,i) = 0; + intra15(i,:) = 0; + %assert(~any(bnet.inter(i,:))) +end +dag15 = [intra15 bnet.inter; + zeros(ss) bnet.intra]; +ns = bnet.node_sizes(:); +ns(nonint) = 1; % disconnected nodes get size 1 +obs_nodes = [bnet.observed(:) bnet.observed(:)+ss]; +bnet15 = mk_bnet(dag15, ns, 'discrete', bnet.dnodes, 'equiv_class', bnet.equiv_class(:), ... + 'observed', obs_nodes(:)); diff --git a/sourcecodes/bnt-master/BNT/general/noisyORtoTable.m b/sourcecodes/bnt-master/BNT/general/noisyORtoTable.m new file mode 100644 index 00000000..18dadd46 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/noisyORtoTable.m @@ -0,0 +1,32 @@ +function CPT = noisyORtoTable(inhibit, leak_inhibit) +% NOISYORTOTABLE Convert noisyOR distribution to CPT +% function CPT = noisyORtoTable(inhibit, leak_inhibit) +% +% inhibit(i) = prob i'th parent will be inhibited (flipped from 1 to 0) +% leak_inhibit - optional suppression of leak +% CPT(U1,...,Un, X) = Pr(X|U1,...,Un) where the Us are the parents (excluding leak). +% State 1 = off, 2 = on + +if nargin < 2, leak_inhibit = 1; end + +q = [leak_inhibit inhibit(:)']; + +if length(q)==1 + CPT = [q 1-q]; + return; +end + +n = length(q); +Bn = ind2subv(2*ones(1,n), 1:(2^n))-1; % all n bit vectors, with the left most column toggling fastest (LSB) +CPT = zeros(2^n, 2); +% Pr(X=0 | U_1 .. U_n) = prod_{i: U_i = on} q_i = prod_i q_i ^ U_i = exp(u' * log(q_i)) +% This method is problematic when q contains zeros + +Q = repmat(q(:)', 2^n, 1); +Q(logical(~Bn)) = 1; +CPT(:,1) = prod(Q,2); +CPT(:,2) = 1-CPT(:,1); + +CPT = reshape(CPT(2:2:end), 2*ones(1,n)); % skip cases in which the leak is off + + diff --git a/sourcecodes/bnt-master/BNT/general/partition_dbn_nodes.m b/sourcecodes/bnt-master/BNT/general/partition_dbn_nodes.m new file mode 100644 index 00000000..5185c6cb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/partition_dbn_nodes.m @@ -0,0 +1,17 @@ +function [pnodes, tnodes] = partition_dbn_nodes(intra, inter) +% PARTITION_DBN_NODES Divide the nodes into a DBN into persistent and transient. +% [pnodes, tnodes] = partition_dbn_nodes(intra, inter) +% Persistent nodes have children in the next time slice, transient nodes do not. + +ss = length(intra); +pnodes = []; +tnodes = []; +for i=1:ss + cs = children(inter, i); + if isempty(cs) + tnodes = [tnodes i]; + else + pnodes = [pnodes i]; + end +end + diff --git a/sourcecodes/bnt-master/BNT/general/partition_matrix_vec_3.m b/sourcecodes/bnt-master/BNT/general/partition_matrix_vec_3.m new file mode 100644 index 00000000..86873a3e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/partition_matrix_vec_3.m @@ -0,0 +1,22 @@ +function [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(A, B, C, n1, n2, bs) + +dom = myunion(n1, n2); +n1i = block(find_equiv_posns(n1, dom), bs(dom)); +n2i = block(find_equiv_posns(n2, dom), bs(dom)); + + + A1 = A(n1i); + A2 = A(n2i); + if isempty(B) + B1 = zeros(size(n1i, 2),size(B, 2)); + B2 = zeros(size(n2i, 2),size(B, 2)); + else + B1 = B(n1i, :); + B2 = B(n2i, :); + end + + + C11 = C(n1i, n1i); + C12 = C(n1i, n2i); + C21 = C(n2i, n1i); + C22 = C(n2i, n2i); diff --git a/sourcecodes/bnt-master/BNT/general/sample_bnet.m b/sourcecodes/bnt-master/BNT/general/sample_bnet.m new file mode 100644 index 00000000..3c219c24 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/sample_bnet.m @@ -0,0 +1,34 @@ +function sample = sample_bnet(bnet, varargin) +% SAMPLE_BNET Generate a random sample from a Bayes net. +% SAMPLE = SAMPLE_BNET(BNET, ...) +% +% sample{i} contains the value of the i'th node. +% i.e., the result is an Nx1 cell array. +% Nodes are sampled in the order given by bnet.order. +% +% Optional arguments: +% +% evidence - initial evidence; if evidence{i} is non-empty, node i won't be sampled. + +% set defauly params +n = length(bnet.dag); +sample = cell(n,1); + +% get optional params +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'evidence', sample = args{i+1}(:); + otherwise, error(['unrecognized argument ' args{i}]) + end +end + +for j=bnet.order(:)' + if isempty(sample{j}) + %ps = parents(bnet.dag, j); + ps = bnet.parents{j}; + e = bnet.equiv_class(j); + sample{j} = sample_node(bnet.CPD{e}, sample(ps)); + end +end diff --git a/sourcecodes/bnt-master/BNT/general/sample_bnet_nocell.m b/sourcecodes/bnt-master/BNT/general/sample_bnet_nocell.m new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/sample_bnet_nocell.m diff --git a/sourcecodes/bnt-master/BNT/general/sample_dbn.m b/sourcecodes/bnt-master/BNT/general/sample_dbn.m new file mode 100644 index 00000000..0477ac4a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/sample_dbn.m @@ -0,0 +1,73 @@ +function seq = sample_dbn(bnet, varargin) +% SAMPLE_DBN Generate a random sequence from a DBN. +% seq = sample_dbn(bnet, ...) +% +% seq{i,t} contains the values of the i'th node in the t'th slice. +% +% Optional arguments: +% +% length - length of sequence to be generated (can also just use sample_dbn(bnet,T)) +% stop_test - name of a function which is used to decide when to stop; +% This will be called as feval(stop_test, seq(:,t)) +% i.e., stop_test is passed a cell array containing all the nodes in the current slice. +% evidence - initial evidence; if evidence{i,t} is non-empty, this node won't be sampled. + +args = varargin; +nargs = length(args); + +if (nargs == 1) & ~isstr(args{1}) + % Old syntax: sample_dbn(bnet, T) + T = args{1}; +else + % get length + T = 1; + for i=1:2:nargs + switch args{i}, + case 'length', T = args{i+1}; + case 'evidence', T = size(args{i+1}, 2); + end + end +end + +ss = length(bnet.intra); +% set default arguments +seq = cell(ss, T); +stop_test = []; +for i=1:2:nargs + switch args{i}, + case 'evidence', seq = args{i+1}; % initialise observed nodes + case 'stop_test', stop_test = args{i+1}; + end +end + +t = 1; +for i=1:ss + if ~isempty(stop_test) | isempty(seq{i,t}) + ps = parents(bnet.dag, i); + e = bnet.equiv_class(i,1); + pvals = seq(ps); + seq{i,t} = sample_node(bnet.CPD{e}, pvals); + %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)' + end +end +t = 2; +done = 0; +while ~done + for i=1:ss + if ~isempty(stop_test) | isempty(seq{i,t}) + ps = parents(bnet.dag, i+ss) + (t-2)*ss; + e = bnet.equiv_class(i,2); + pvals = seq(ps); + seq{i,t} = sample_node(bnet.CPD{e}, pvals); + %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)' + end + end + if ~isempty(stop_test) + done = feval(stop_test, seq(:,t)); + else + if t==T + done = 1; + end + end + t = t + 1; +end diff --git a/sourcecodes/bnt-master/BNT/general/score_bnet_complete.m b/sourcecodes/bnt-master/BNT/general/score_bnet_complete.m new file mode 100644 index 00000000..93caf1d7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/score_bnet_complete.m @@ -0,0 +1,28 @@ +function L = log_lik_complete(bnet, cases, clamped) +% LOG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m), theta_i) for a completely observed data set +% L = log_lik_complete(bnet, cases, clamped) +% +% If there is a missing data, you must use an inference engine. +% cases(i,m) is the value assigned to node i in case m. +% (If there are vector-valued nodes, cases should be a cell array.) +% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros) +% Clamped nodes contribute a factor of 1.0 to the likelihood. + +if iscell(cases), usecell = 1; else usecell = 0; end + +n = length(bnet.dag); +ncases = size(cases, 2); +if n ~= size(cases, 1) + error('data should be of size nnodes * ncases'); +end + +if nargin < 3, clamped = zeros(n,ncases); end + +L = 0; +for i=1:n + ps = parents(bnet.dag, i); + e = bnet.equiv_class(i); + u = find(clamped(i,:)==0); + L = L + log_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u)); +end + diff --git a/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_gaussian.m b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_gaussian.m new file mode 100644 index 00000000..8d4a3481 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_gaussian.m @@ -0,0 +1,12 @@ +function marg2 = shrink_obs_dims_in_gaussian(marg1, dom, evidence, ns) +% SHRINK_OBS_DIMS_IN_GAUSSIAN Remove observed dimensions from mu/Sigma +% function marg2 = shrink_obs_dims_in_gaussian(marg1, dom, evidence, ns) + +% This is used by loopy + +hdom = dom(isemptycell(evidence(dom))); +ndx = find_equiv_posns(hdom, dom); +b = block(ndx, ns(dom)); +marg2.mu = marg1.mu(b); +marg2.Sigma = marg1.Sigma(b,b); +marg2.domain = marg1.domain; diff --git a/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_table.m b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_table.m new file mode 100644 index 00000000..37c7d804 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_table.m @@ -0,0 +1,15 @@ +function T2 = shrink_obs_dims_in_table(T1, dom, evidence) +% SHRINK_OBS_DIMS_IN_TABLE Set observed dimensions to size 1 +% T2 = shrink_obs_dims_in_table(T1, dom, evidence) +% +% If 'T1' contains observed nodes, it will have 0s in the positions that are +% inconsistent with the evidence. We now remove these 0s and set the corresponding dimensions to +% size 1, to be consistent with the way most inference engines handle evidence, which is to +% shrink observed nodes before doing inference. + +% This is used by pearl and enumerative inf. engines. + +odom = dom(~isemptycell(evidence(dom))); +vals = cat(1,evidence{odom}); +ndx = mk_multi_index(length(dom), find_equiv_posns(odom, dom), vals(:)); +T2 = T1(ndx{:}); diff --git a/sourcecodes/bnt-master/BNT/general/solve_limid.m b/sourcecodes/bnt-master/BNT/general/solve_limid.m new file mode 100644 index 00000000..0009cd7c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/solve_limid.m @@ -0,0 +1,68 @@ +function [strategy, MEU, niter] = solve_limid(engine, varargin) +% SOLVE_LIMID Find the (locally) optimal strategy for a LIMID +% [strategy, MEU, niter] = solve_limid(inf_engine, ...) +% +% strategy{d} = stochastic policy for node d (a decision node) +% MEU = maximum expected utility +% niter = num iterations used +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default in brackets] +% +% max_iter - max. num. iterations [ 1 ] +% tol - tolerance required of consecutive MEU values, used to assess convergence [1e-3] +% order - order in which decision nodes are optimized [ reverse numerical order ] +% +% e.g., solve_limid(engine, 'tol', 1e-2, 'max_iter', 10) + +bnet = bnet_from_engine(engine); + +% default values +max_iter = 1; +tol = 1e-3; +D = bnet.decision_nodes; +order = D(end:-1:1); + +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 'order', order = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end +end + +CPDs = bnet.CPD; +ns = bnet.node_sizes; +N = length(ns); +evidence = cell(1,N); +strategy = cell(1, N); + +iter = 1; +converged = 0; +oldMEU = 0; +while ~converged & (iter <= max_iter) + for d=order(:)' + engine = enter_evidence(engine, evidence, 'exclude', d); + [m, pot] = marginal_family(engine, d); + %pot = marginal_family_pot(engine, d); + [policy, score] = upot_to_opt_policy(pot); + e = bnet.equiv_class(d); + CPDs{e} = set_fields(CPDs{e}, 'policy', policy); + engine = update_engine(engine, CPDs); + strategy{d} = policy; + end + engine = enter_evidence(engine, evidence); + [m, pot] = marginal_nodes(engine, []); + %pot = marginal_family_pot(engine, []); + [dummy, MEU] = upot_to_opt_policy(pot); + if approxeq(MEU, oldMEU, tol) + converged = 1; + end + oldMEU = MEU; + iter = iter + 1; +end +niter = iter - 1; diff --git a/sourcecodes/bnt-master/BNT/general/unroll_dbn_topology.m b/sourcecodes/bnt-master/BNT/general/unroll_dbn_topology.m new file mode 100644 index 00000000..6b9ddcdb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/unroll_dbn_topology.m @@ -0,0 +1,28 @@ +function M = unroll_dbn_topology(intra, inter, T, intra1) +% UNROLL_DBN_TOPOLOGY Make the block diagonal adjacency matrix for a DBN consisting of T slices +% M = unroll_dbn_topology(intra, inter, T, intra1) +% +% intra is the connectivity within a slice, inter between two slices. +% M will have intra along the diagonal, and inter one above the diagonal. +% intra1 is an optional argumnet, in case the intra is different for the first slice. + +if nargin < 4, intra1 = intra; end + +ss = length(intra); % slice size +M = sparse(ss*T, ss*T); + +b = 1:ss; +M(b,b) = intra1; +M(b,b+ss) = inter; + +for t=2:T-1 + b = (1:ss) + (t-1)*ss; + M(b,b) = intra; + M(b,b+ss) = inter; +end + +t = T; +b = (1:ss) + (t-1)*ss; +M(b,b) = intra; + + diff --git a/sourcecodes/bnt-master/BNT/general/unroll_higher_order_topology.m b/sourcecodes/bnt-master/BNT/general/unroll_higher_order_topology.m new file mode 100644 index 00000000..29e78ef2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/unroll_higher_order_topology.m @@ -0,0 +1,30 @@ +function M = unroll_higher_order_topology(intra, inter, T, intra1) +% UNROLL_DBN_TOPOLOGY Make the block diagonal adjacency matrix for a DBN consisting of T slices +% M = unroll_dbn_topology(intra, inter, T, intra1) +% +% intra is the connectivity within a slice, inter between two slices. +% M will have intra along the diagonal, and inter one above the diagonal. +% intra1 is an optional argumnet, in case the intra is different for the first slice. + +if nargin < 4 + intra1 = intra; +end; + + +ss = length(intra); % slice size +M = sparse(ss*T, ss*T); +[rows,columns,order] = size(inter); +for t1 = 1:T + b = 1 + (t1 - 1)*ss : t1*ss; + if t1 == 1 + M(b,b) = intra1; + else + M(b,b) = intra; + end + for t2 = 1:order + if t1 + t2 <= T + M(b,b+t2*ss) = inter(:,:,t2); + end + end +end + diff --git a/sourcecodes/bnt-master/BNT/general/unroll_set.m b/sourcecodes/bnt-master/BNT/general/unroll_set.m new file mode 100644 index 00000000..277efcd2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/general/unroll_set.m @@ -0,0 +1,7 @@ +function U = unroll_set(S, ss, T) +% UNROLL_SET Make T shifted copies of the set of nodes S in a slice of size ss. +% U = unroll_set(S, ss, T) + +offset = repmat(0:ss:(T-1)*ss, [length(S) 1]); +U = repmat(S(:), [1 T]) + offset; + |