diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine')
21 files changed, 764 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Entries new file mode 100644 index 00000000..c19ebdd4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Entries @@ -0,0 +1,4 @@ +/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/gibbs_sampling_inf_engine.m/1.1.1.1/Wed May 29 15:59:56 2002// +/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002// +D/private//// diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Repository new file mode 100644 index 00000000..3338daf9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/static/@gibbs_sampling_inf_engine diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/enter_evidence.m new file mode 100644 index 00000000..0710d5c8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/enter_evidence.m @@ -0,0 +1,29 @@ +function [engine, loglik] = enter_evidence(engine, evidence) +% ENTER_EVIDENCE Add the specified evidence to the network (gibbs_sampling_inf_engine) +% [engine, loglik] = enter_evidence(engine, evidence) +% +% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value +% +% loglik is not computed... we just return a 0 value + +bnet = bnet_from_engine(engine); + +engine.hnodes = find(isemptycell(evidence)); +engine.onodes = mysetdiff(1:length(evidence), engine.hnodes); + +engine.evidence = zeros(engine.slice_size, 1); + +% Reset all counts since they are no longer valid +engine.marginal_counts = {}; +%engine.state = sample_bnet (bnet, 1, 0); +engine.state = cell2num(sample_bnet(bnet)); + +% For speed, we use a normal (not cell) array. We're making use of +% the current restriction to discrete nodes. +for i = engine.onodes + engine.evidence(i) = evidence{i}; +end + +loglik = 0; + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/gibbs_sampling_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/gibbs_sampling_inf_engine.m new file mode 100644 index 00000000..3dc4b361 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/gibbs_sampling_inf_engine.m @@ -0,0 +1,104 @@ +function engine = gibbs_sampling_inf_engine(bnet, varargin) +% GIBBS_SAMPLING_INF_ENGINE +% +% engine = gibbs_sampling_inf_engine(bnet, ...) +% +% Optional parameters [default in brackets] +% 'burnin' - How long before you start using the samples [100]. +% 'gap' - how often you use the samples in the estimate [1]. +% 'T' - number of samples [1000] +% i.e, number of node flips (so, for +% example if there are 10 nodes in the bnet, and T is 1000, each +% node will get flipped 100 times (assuming a deterministic schedule)) +% The total running time is proportional to burnin + T*gap. +% +% 'order' - if the sampling schedule is deterministic, use this +% parameter to specify the order in which nodes are sampled. +% Order is allowed to include multiple copies of nodes, which is +% useful if you want to, say, focus sampling on particular nodes. +% Default is to use a deterministic schedule that goes through the +% nodes in order. +% +% 'sampling_dist' - when using a stochastic sampling method, at +% each step the node to sample is chosen according to this +% distribution (may be unnormalized) +% +% The sampling_dist and order parameters shouldn't both be used, +% and this will cause an assert. +% +% +% Written by "Bhaskara Marthi" <bhaskara@cs.berkeley.edu> Feb 02. + + +engine.burnin = 100; +engine.gap = 1; +engine.T = 1000; +use_default_order = 1; +engine.deterministic = 1; +engine.order = {}; +engine.sampling_dist = {}; + +if nargin >= 2 + args = varargin; + nargs = length(args); + for i = 1:2:nargs + switch args{i} + case 'burnin' + engine.burnin = args{i+1}; + case 'gap' + engine.gap = args{i+1}; + case 'T' + engine.T = args{i+1}; + case 'order' + assert (use_default_order); + use_default_order = 0; + engine.order = args{i+1}; + case 'sampling_dist' + assert (use_default_order); + use_default_order = 0; + engine.deterministic = 0; + engine.sampling_dist = args{i+1}; + otherwise + error(['unrecognized parameter to gibbs_sampling_inf_engine']); + end + end +end + +engine.slice_size = size(bnet.dag, 2); +if (use_default_order) + engine.order = 1:engine.slice_size; +end +engine.hnodes = []; +engine.onodes = []; +engine.evidence = []; +engine.state = []; +engine.marginal_counts = {}; + +% Precompute the strides for each CPT +engine.strides = compute_strides(bnet); + +% Precompute graphical information +engine.families = compute_families(bnet); +engine.children = compute_children(bnet); + +% For convenience, store the CPTs as tables rather than objects +engine.CPT = get_cpts(bnet); + +engine = class(engine, 'gibbs_sampling_inf_engine', inf_engine(bnet)); + + + + + + + + + + + + + + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/marginal_nodes.m new file mode 100644 index 00000000..8df75552 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/marginal_nodes.m @@ -0,0 +1,135 @@ +function [marginal, engine] = marginal_nodes(engine, nodes, varargin); +% MARGINAL_NODES Compute the marginal on the specified query nodes +% (gibbs_sampling_engine) +% [marginal, engine] = marginal_nodes(engine, nodes, ...) +% +% returns Pr(X(nodes) | X(observedNodes)) +% +% The engine is also modified, and so it is returned as well, since +% Matlab doesn't support passing by reference(!) So +% if you want to, for example, incrementally run gibbs for a few 100 +% steps at a time, you should use the returned value. +% +% Optional arguments : +% +% 'reset_counts' is 1 if you want to reset the counts made in the +% past, and 0 otherwise (if the current query nodes are different +% from the previous query nodes, or if marginal_nodes has not been +% called before, reset_counts should be set to 1). +% By default it is 1. + + +reset_counts = 1; + +if (nargin > 3) + args = varargin; + nargs = length(args); + for i = 1:2:nargs + switch args{i} + case 'reset_counts' + reset_counts = args{i+1}; + otherwise + error(['Incorrect argument to gibbs_sampling_engine/' ... + ' marginal_nodes']); + end + end +end + +% initialization stuff +bnet = bnet_from_engine(engine); +slice_size = engine.slice_size; +hnodes = engine.hnodes; +onodes = engine.onodes; +nonqnodes = mysetdiff(1:slice_size, nodes); +gap = engine.gap; +burnin = engine.burnin; +T_max = engine.T; +ns = bnet.node_sizes(nodes); + + +% Cache the strides for the marginal table +marg_strides = [1 cumprod(ns(1:end-1))]; + +% Reset counts if necessary +if (reset_counts == 1) + %state = sample_bnet(bnet, 1, 0); + %state = cell2num(sample_bnet(bnet, 'evidence', num2cell(engine.evidence))); + state = cell2num(sample_bnet(bnet)); + state(onodes) = engine.evidence(onodes); + if (length(ns) == 1) + marginal_counts = zeros(ns(1),1); + else + marginal_counts = zeros(ns); + end + +% Otherwise, use the counts that have been stored in the engine +else + state = engine.state; + state(onodes, :) = engine.evidence(onodes, :); + marginal_counts = engine.marginal_counts; +end + +if (engine.deterministic == 1) + pos = 1; + order = engine.order; + orderSize = length(engine.order); +else + sampling_dist = normalise(engine.sampling_dist); +end + + +for t = 1:(T_max*gap+burnin) + + % First, select node m to sample + if (engine.deterministic == 1) + m = engine.order(pos); + pos = pos+1; + if (pos > orderSize) + pos = 1; + end + else + m = my_sample_discrete(sampling_dist); + end + + + % If the node is observed, then don't bother resampling + if (myismember(m, onodes)) + continue; + end + + % Next, compute the posterior + post = compute_posterior (bnet, state, m, engine.strides, engine.families, ... + engine.children, engine.CPT); + state(m) = my_sample_discrete(post); + + % Now update our monte carlo estimate of the posterior + % distribution on the query node + if ((mod(t-burnin, gap) == 0) & (t > burnin)) + + vals = state(nodes); + index = 1+marg_strides*(vals-1); + marginal_counts(index) = marginal_counts(index)+1; + end +end + +% Store results for future computation. Note that we store +% unnormalized counts +engine.state = state; +engine.marginal_counts = marginal_counts; + +marginal.T = normalise(marginal_counts); + + + + + + + + + + + + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m new file mode 100644 index 00000000..772f137c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m @@ -0,0 +1,5 @@ +function c = CPT(bnet, i) +% CPT Helper function avoid having to type in +% CPD_to_CPT(bnet.CPD{i}) every time + +c = CPD_to_CPT(bnet.CPD{i}); \ No newline at end of file diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries new file mode 100644 index 00000000..0919a694 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries @@ -0,0 +1,13 @@ +/CPT.m/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_children.m/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_families.m/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_families_dbn.m/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_posterior.c/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_posterior_dbn.m/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_strides.m/1.1.1.1/Wed May 29 15:59:56 2002// +/get_cpts.m/1.1.1.1/Wed May 29 15:59:56 2002// +/get_slice_dbn.c/1.1.1.1/Wed May 29 15:59:56 2002// +/get_slice_dbn.m/1.1.1.1/Wed May 29 15:59:56 2002// +/my_sample_discrete.m/1.1.1.1/Wed May 29 15:59:56 2002// +/sample_single_discrete.c/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository new file mode 100644 index 00000000..a3027631 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/static/@gibbs_sampling_inf_engine/private diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_children.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_children.m new file mode 100644 index 00000000..3af799f8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_children.m @@ -0,0 +1,12 @@ +function c = compute_children(bnet) +% COMPUTE_CHILDREN +% precomputes the children of nodes in a bnet +% +% The return value is a cell array for now + +ss = size(bnet.dag, 1); +c = cell(ss, 1); +for i = 1:ss + c{i} = children(bnet.dag, i); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m new file mode 100644 index 00000000..e75974cc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m @@ -0,0 +1,12 @@ +function families = compute_families(bnet) +% COMPUTE_FAMILIES +% precomputes the families of nodes in a bnet +% +% The return value is a cell array for now + +ss = size(bnet.dag, 1); +families = cell(ss, 1); +for i = 1:ss + families{i} = family(bnet.dag, i); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m new file mode 100644 index 00000000..7647bc28 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m @@ -0,0 +1,13 @@ +function families = compute_families_dbn(bnet) +% COMPUTE_FAMILIES +% precomputes the families of nodes in a dbn +% +% The return value is a cell array for now + +ss = size(bnet.intra, 1); +families = cell(ss, 2); +for i = 1:ss + families{i, 1} = family(bnet.dag, i, 1); + families{i, 2} = family(bnet.dag, i, 2); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c new file mode 100644 index 00000000..3c61b7f3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c @@ -0,0 +1,107 @@ +#include "mex.h" + +/* Helper function that extracts a one-dimensional slice from a cpt */ +/* +void multiplySlice(mxArray *bnet, mxArray *state, int i, int nsi, int j, + mxArray *strides, mxArray *fam, mxArray *cpts, + double *y) +*/ +void multiplySlice(const mxArray *bnet, const mxArray *state, int i, int nsi, int j, + const mxArray *strides, const mxArray *fam, const mxArray *cpts, + double *y) +{ + mxArray *ec, *cpt, *family; + double *ecElts, *cptElts, *famElts, *strideElts, *ev; + int c1, k, famSize, startInd, strideStride, pos, stride; + + strideStride = mxGetM(strides); + strideElts = mxGetPr(strides); + + ev = mxGetPr(state); + + /* Get the CPT */ + ec = mxGetField (bnet, 0, "equiv_class"); + ecElts = mxGetPr(ec); + k = (int) ecElts[j-1]; + cpt = mxGetCell (cpts, k-1); + cptElts = mxGetPr (cpt); + + /* Get the family vector for this cpt */ + family = mxGetCell (fam, j-1); + famSize = mxGetNumberOfElements (family); + famElts = mxGetPr (family); + + /* Figure out starting position and stride */ + startInd = 0; + for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos +=strideStride) { + if (famElts[c1] != i) { + startInd += strideElts[pos]*(ev[(int)famElts[c1]-1]-1); + } + else { + stride = strideElts[pos]; + } + } + + for (c1 = 0, pos = startInd; c1 < nsi; c1++, pos+=stride) { + y[c1] *= cptElts[pos]; + } +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray + *prhs[]) +{ + double *pi, *nsElts, *y, *childrenElts; + mxArray *ns, *children; + double sum; + int i, nsi, c1, numChildren; + + pi = mxGetPr(prhs[2]); + i = (int) pi[0]; + + ns = mxGetField(prhs[0], 0, "node_sizes"); + nsElts = mxGetPr(ns); + nsi = (int) nsElts[i-1]; + + /* Initialize the posterior */ + plhs[0] = mxCreateDoubleMatrix (1, nsi, mxREAL); + y = mxGetPr(plhs[0]); + for (c1 = 0; c1 < nsi; c1++) { + y[c1] = 1; + } + + /* Multiply in the cpt of the node i */ + multiplySlice(prhs[0], prhs[1], i, nsi, i, prhs[3], prhs[4], + prhs[6], y); + + + /* Multiply in cpts of children of i */ + children = mxGetCell (prhs[5], i-1); + numChildren = mxGetNumberOfElements (children); + childrenElts = mxGetPr (children); + + for (c1 = 0; c1 < numChildren; c1++) { + int j; + j = (int) childrenElts[c1]; + multiplySlice (prhs[0], prhs[1], i, nsi, j, prhs[3], prhs[4], + prhs[6], y); + } + + sum = 0; + /* normalize! */ + for (c1 = 0; c1 < nsi; c1++) { + sum += y[c1]; + } + + for (c1 = 0; c1 < nsi; c1++) { + y[c1] /= sum; + } +} + + + + + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m new file mode 100644 index 00000000..e9a69b24 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m @@ -0,0 +1,59 @@ +function post = compute_posterior_dbn(bnet, state, i, n, strides, families, ... + CPT) +% COMPUTE_POSTERIOR +% +% post = compute_posterior(bnet, state, i, n, strides, families, +% cpts) +% +% Compute the posterior distribution on node X_i^n of a DBN, +% conditional on evidence in the cell array state +% +% strides is the cached result of compute_strides(bnet) +% families is the cached result of compute_families(bnet) +% cpt is the cached result of get_cpts(bnet) +% +% post is a one-dimensional table + + + +% First multiply in the cpt of the node itself +post = get_slice_dbn(bnet, state, i, n, i, n, strides, families, CPT); +post = post(:); + +% Then multiply in CPTs of children that are in this slice +for j = children(bnet.intra, i) + slice = get_slice_dbn(bnet, state, j, n, i, n, strides, families, CPT); + post = post.*slice(:); +end + +% Finally, if necessary, multiply in CPTs of children in the next +% slice +if (n < size(state,2)) + for j = children(bnet.inter, i) + slice = get_slice_dbn(bnet, state, j, n+1, i, n, strides, families, ... + CPT); + post = post.*slice(:); + end +end + +post = normalise(post); + + + + + + + + + + + + + + + + + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m new file mode 100644 index 00000000..a8e26c25 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m @@ -0,0 +1,27 @@ +function strides = compute_strides(bnet) +% COMPUTE_STRIDES For each CPT and each variable in that CPT, +% returns the stride of that variable. So in future, we can +% quickly extract a slice of the CPT. +% +% The return value is a 2d array, where strides(i,j) contains the +% stride of the jth variable in the ith CPT. Cell arrays would +% have saved space but they are slower. +% + +num_cpts = size(bnet.CPD, 2); +max_cpt_dim = 1 + max(sum(bnet.dag)); +strides = zeros(num_cpts, max_cpt_dim); + +for i = 1:num_cpts + c = CPT(bnet, i); + siz = size(CPT(bnet, i)); + + % Deal with the special case of a 1-d array separately + if siz(2) == 1 + dim = 1; + else + dim = size(siz, 2); + end + + strides(i, 1:dim ) = [1 cumprod(siz(1:dim-1))]; +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m new file mode 100644 index 00000000..77c86070 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m @@ -0,0 +1,8 @@ +function c = get_cpts(bnet) +% Get all the cpts in tabular form + +cpds = bnet.CPD; +c = cell(size(cpds)); +for i = 1:length(c) + c{i} = CPT(bnet, i); +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c new file mode 100644 index 00000000..33540eff --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c @@ -0,0 +1,116 @@ +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray + *prhs[]) +{ + double *pn, *pi, *pj, *pm, *y, *ecElts, *pcpt, *famElts, *strideElts, + *ev, *nsElts; + int i, k, j, m, n; + mxArray *ec, *cpt, *fam, *ns; + int c1, famSize, nsj; + int strideStride, startInd, stride, pos, numNodes; + + const int BNET = 0; + const int STATE = 1; + const int STRIDES = 6; + const int FAMILIES = 7; + const int CPT = 8; + + pn = mxGetPr(prhs[3]); + n = (int) pn[0]; + pi = mxGetPr(prhs[2]); + i = (int) pi[0]; + pj = mxGetPr(prhs[4]); + j = (int) pj[0]; + pm = mxGetPr(prhs[5]); + m = (int) pm[0]; + ev = mxGetPr(prhs[STATE]); + ns = mxGetField (prhs[BNET], 0, "node_sizes"); + nsElts = mxGetPr (ns); + numNodes = mxGetM(ns); + + strideStride = mxGetM(prhs[STRIDES]); + strideElts = mxGetPr(prhs[STRIDES]); + + + + /* Treat the case n = 1 separately */ + if (pn[0] == 1) { + + /* Get the appropriate CPT */ + ec = mxGetField (prhs[BNET], 0, "eclass1"); + ecElts = mxGetPr(ec); + k = (int) ecElts[i-1]; + cpt = mxGetCell (prhs[8], k-1); + pcpt = mxGetPr(cpt); + + nsj = (int) nsElts[j-1]; + + /* Get the correct family vector */ + /* (Note : MEX is painful) */ + fam = mxGetCell (prhs[FAMILIES], i - 1); + famSize = mxGetNumberOfElements(fam); + famElts = mxGetPr(fam); + + + /* Figure out starting position and stride */ + startInd = 0; + for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos+=strideStride) { + if (famElts[c1] != j) { + startInd += strideElts[pos]*(ev[(int)famElts[c1]-1]-1); + } + else { + stride = strideElts[pos]; + } + } + + plhs[0] = mxCreateDoubleMatrix (1, nsj, mxREAL); + y = mxGetPr(plhs[0]); + for (c1 = 0, pos = startInd; c1 < nsj; c1++, pos+=stride) { + y[c1] = pcpt[pos]; + } + } + + /* Handle the case n > 1 */ + else { + + /* Get the appropriate CPT */ + ec = mxGetField (prhs[BNET], 0, "eclass2"); + ecElts = mxGetPr(ec); + k = (int) ecElts[i-1]; + cpt = mxGetCell (prhs[8], k-1); + pcpt = mxGetPr(cpt); + + /* Figure out size of slice */ + if (m == 1) { + nsj = (int) nsElts[j-1]; + } + else { + nsj = (int) nsElts[j-1+numNodes]; + } + + /* Figure out family */ + fam = mxGetCell (prhs[FAMILIES], i - 1 + numNodes); + famSize = mxGetNumberOfElements(fam); + famElts = mxGetPr(fam); + + startInd = 0; + for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos+=strideStride) { + int f = (int) famElts[c1]; + + if (((f == j+numNodes) && (m == n)) || ((f == j) && (m == + n-1))) { + stride = strideElts[pos]; + } + else { + startInd += strideElts[pos] * (ev[f-1+((n-2)*numNodes)]-1); + } + } + + plhs[0] = mxCreateDoubleMatrix(1,nsj, mxREAL); + y = mxGetPr(plhs[0]); + for (c1 = 0, pos = startInd; c1 < nsj; c1++, pos+=stride) { + y[c1] = pcpt[pos]; + } + } +} diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m new file mode 100644 index 00000000..22841784 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m @@ -0,0 +1,87 @@ +function slice = get_slice_dbn(bnet, state, i, n, j, m, strides, families, ... + CPT) +% slice = get_slice(bnet, state, i, n, j, m, strides, families, cpt) +% +% GET_SLICE get one-dimensional slice of the CPT for node X_i^n +% that corresponds to the different values of X_j^m, where all +% other nodes have values given by state. +% strides is the result of +% calling compute_strides(bnet) +% families is the result of calling compute_families(bnet) +% cpts is the result of calling get_cpts(bnet) +% +% slice is a 1-d array + + +if (n == 1) + + k = bnet.eclass1(i); + c = CPT{k}; + + % Figure out evidence on family + fam = families{i, 1}; + ev = state(fam, 1); + + % Remove evidence on node j + pos = find(fam == j); + ev(pos) = 1; + dim = size(ev, 1); + + % Compute initial index and stride + start_ind = 1+strides(k, 1:dim)*(ev-1); + stride = strides(k, pos); + + % Compute the slice + slice = c(start_ind:stride:start_ind+(bnet.node_sizes(j, 1)-1)*stride); + +else + + k = bnet.eclass2(i); + c = CPT{k}; + + fam = families{i, 2}; + ss = length(bnet.intra); + + % Divide the family into nodes in this time step and nodes in the + % previous time step + this_time_step = fam(find(fam > ss)); + prev_time_step = fam(find(fam <= ss)); + + % Normalize the node numbers + this_time_step = this_time_step - ss; + + % Get the evidence + this_step_ev = state(this_time_step, n); + prev_step_ev = state(prev_time_step, n-1); + + % Remove the evidence for X_j^m + if (m == n) + pos = find(this_time_step == j); + this_step_ev(pos) = 1; + pos = pos + size(prev_time_step, 2); + else + assert (m == n-1); + pos = find(prev_time_step == j); + prev_step_ev(pos) = 1; + end + + % Combine the two time steps + ev = [prev_step_ev; this_step_ev]; + dim = size(ev, 1); + + + % Compute starting index and stride + start_ind = 1 + strides(k, 1:dim)*(ev-1); + stride = strides(k, pos); + + % Compute slice + if (m == 1) + q = 1; + else + q = 2; + end + slice = c(start_ind:stride:start_ind+(bnet.node_sizes(j, q)-1)*stride); +end + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m new file mode 100644 index 00000000..70f0615b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m @@ -0,0 +1,7 @@ +function M = my_sample_discrete(prob) +% A faster version that calls a c subfunction. Will update one +% day to have r and c parameters as well + +R = rand (1,1); +M = sample_single_discrete(R, prob); + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c new file mode 100644 index 00000000..36112de6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c @@ -0,0 +1,22 @@ +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray + *prhs[]) +{ + double *y, *pr, *dist; + int k, distSize; + double r, cumSum; + + plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL); + y = mxGetPr (plhs[0]); + + pr = mxGetPr (prhs[0]); + r = pr[0]; + + dist = mxGetPr (prhs[1]); + distSize = mxGetNumberOfElements (prhs[1]); + + for (k = 0, cumSum = 0; (k < distSize) && (r >= cumSum); cumSum += dist[k], k++); + + y[0] = k; +} |