diff options
| author | ziejd2 | 2018-03-14 23:23:33 -0500 |
|---|---|---|
| committer | GitHub | 2018-03-14 23:23:33 -0500 |
| commit | 1ff6baa44e22b91eefb48aea6f3befa078c0489b (patch) | |
| tree | e0fd79d2e32fd2aedda2eadaed0f19af3514c520 /sourcecodes/bnt-master/BNT/potentials | |
| parent | 6882395afdadf4e982b25b5215071a0932730950 (diff) | |
| parent | c80226899f5cdd9f11c163817d59445213f5bef0 (diff) | |
| download | BNW-1ff6baa44e22b91eefb48aea6f3befa078c0489b.tar.gz | |
Merge pull request #1 from ziejd2/octave_php_separate
Octave php separate
Diffstat (limited to 'sourcecodes/bnt-master/BNT/potentials')
145 files changed, 4449 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries new file mode 100644 index 00000000..25f4b469 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries @@ -0,0 +1,15 @@ +/cg_can_to_mom.m/1.1.1.1/Wed May 29 15:59:58 2002// +/cg_mom_to_can.m/1.1.1.1/Wed May 29 15:59:58 2002// +/cgpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/display.m/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/enter_cts_evidence_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/enter_discrete_evidence_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/multiply_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/multiply_pots.m/1.1.1.1/Wed May 29 15:59:58 2002// +/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002// +/set_domain_pot.m/1.1.1.1/Wed Jul 30 13:38:24 2003// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries.Log new file mode 100644 index 00000000..24f16336 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries.Log @@ -0,0 +1 @@ +A D/Old//// diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Repository new file mode 100644 index 00000000..ee9924e1 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@cgpot diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Entries new file mode 100644 index 00000000..680a9046 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Entries @@ -0,0 +1,4 @@ +/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/set_domain_pot.m/1.1.1.1/Wed Jul 30 13:38:08 2003// +/simple_marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Repository new file mode 100644 index 00000000..7feb051e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@cgpot/Old diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/normalize_pot.m new file mode 100644 index 00000000..1debd5bc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/normalize_pot.m @@ -0,0 +1,24 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the CG potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +% Marginalize down to [], so that the normalizing constant becomes Pr(E) +temp = marginalize_pot(cg_can_to_mom(pot), []); +%loglik = temp.mom{1}.logp; +[temp2, loglik] = normalize_pot(temp.mom{1}); + +% Adjust scale factor to reflect the fact that the pot now represents Pr(X | E) instead of Pr(X,E). + +scale = -loglik; +if 1 +switch pot.subtype + case 'm' + for i=1:pot.dsize + pot.mom{i} = rescale_pot(pot.mom{i}, scale); + end + case 'c' + for i=1:pot.dsize + pot.can{i} = rescale_pot(pot.can{i}, scale); + end +end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/set_domain_pot.m new file mode 100644 index 00000000..f0df6ea7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/set_domain_pot.m @@ -0,0 +1,9 @@ +function pot = set_domain_pot(pot, domain) +% SET_DOMAIN_POT Change the domain of a potential (cgpot) +% pot = set_domain_pot(pot, domain) + +delta = domain(1) - pot.domain(1); +assert(all(domain == pot.domain + delta)); +pot.domain = pot.domain + delta; +pot.ddom = pot.ddom + delta; +pot.cdom = pot.cdom + delta; diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/simple_marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/simple_marginalize_pot.m new file mode 100644 index 00000000..a2774636 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/simple_marginalize_pot.m @@ -0,0 +1,86 @@ +function smallpot = marginalize_pot(bigpot, keep) +% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, keep) + +sumover = mysetdiff(bigpot.domain, keep); +csumover = myintersect(sumover, bigpot.cdom); +dsumover = myintersect(sumover, bigpot.ddom); +dkeep = myintersect(keep, bigpot.ddom); +ckeep = myintersect(keep, bigpot.cdom); +%ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer +ns = zeros(1, max(bigpot.domain)); +ns(bigpot.ddom) = bigpot.dsizes; +ns(bigpot.cdom) = bigpot.csizes; + +% sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes. +% Similarly, prod(ns(dsumover))==1 is like isempty(dsumover) + +% Marginalize the cts parts. +% If we are in canonical form, we stay that way, since moment form might not exist. +% Besides, we would like to minimize the number of conversions. +if sum(ns(csumover)) > 0 + if bigpot.subtype == 'm' + for i=1:bigpot.dsize + bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep); + end + else + for i=1:bigpot.dsize + bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep); + end + end +end + +% If we are not marginalizing over any discrete nodes, we are done. +if prod(ns(dsumover))==1 + smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype); + return; +end + +% To marginalize the discrete parts, we must be in moment form. +bigpot = cg_can_to_mom(bigpot); + +I = prod(ns(dkeep)); +J = prod(ns(dsumover)); +C = sum(ns(ckeep)); + +% Reshape bigpot into the form mu1(:,j,i), where i is in dkeep, j is in dsumover +T1 = zeros(I,J); +mu1 = zeros(C,J,I); +Sigma1 = zeros(C,C,J,I); +sum_map = find_equiv_posns(dsumover, bigpot.ddom); +keep_map = find_equiv_posns(dkeep, bigpot.ddom); +iv = zeros(1, length(bigpot.ddom)); % index vector +for i=1:I + keep_iv = ind2subv(ns(dkeep), i); + iv(keep_map) = keep_iv; + for j=1:J + sum_iv = ind2subv(ns(dsumover), j); + iv(sum_map) = sum_iv; + k = subv2ind(ns(bigpot.ddom), iv); + mom = struct(bigpot.mom{k}); % violate object privacy + T1(i,j) = exp(mom.logp); + if C > 0 % so mu1 and Sigma1 are non-empty + mu1(:,j,i) = mom.mu; + Sigma1(:,:,j,i) = mom.Sigma; + end + end +end + +% Collapse the mixture of Gaussians +coef = mk_stochastic(T1); % coef must be convex combination +T2 = sum(T1,2); +T2 = T2 + (T2==0)*eps; +%if C > 0, disp('collapsing onto '); disp(leep); end +mu = []; +Sigma = []; +mom = cell(1,I); +for i=1:I + if C > 0 + [mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:)); + end + logp = log(T2(i)); + mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma); +end + +smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm'); + diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_can_to_mom.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_can_to_mom.m new file mode 100644 index 00000000..951dc910 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_can_to_mom.m @@ -0,0 +1,10 @@ +function pot = cg_can_to_mom(pot) +% CG_CAN_TO_MOM Convert a CG potential from canonical to moment form, if necessary. +% pot = cg_can_to_mom(pot) + +if pot.subtype ~= 'm' + for i=1:pot.dsize + pot.mom{i} = cpot_to_mpot(pot.can{i}); + end + pot.subtype = 'm'; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_mom_to_can.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_mom_to_can.m new file mode 100644 index 00000000..61aa053f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_mom_to_can.m @@ -0,0 +1,10 @@ +function pot = cg_mom_to_can(pot) +% CG_MOM_TO_CAN Convert a CG potential from moment to canonical form, if necessary. +% pot = cg_mom_to_can(pot) + +if pot.subtype ~= 'c' + for i=1:pot.dsize + pot.can{i} = mpot_to_cpot(pot.mom{i}); + end + pot.subtype = 'c'; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/cgpot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cgpot.m new file mode 100644 index 00000000..676e3ce1 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cgpot.m @@ -0,0 +1,38 @@ +function pot = cgpot(ddom, cdom, node_sizes, can, mom, subtype) +% CPOT Make a canonical CG potential. +% function pot = cgpot(ddom, cdom, node_sizes, can, mom, subtype) +% +% node_sizes(i) is the size of the i'th node. +% can and mom default to 0s. +% subtype defaults to 'c'. + +if nargin < 6, subtype = 'c'; end + +pot.ddom = ddom; +pot.cdom = cdom; +node_sizes = node_sizes(:)'; % row vectors print better +pot.domain = myunion(ddom, cdom); +pot.dsizes = node_sizes(pot.ddom); +pot.dsize = prod(node_sizes(pot.ddom)); +pot.csizes = node_sizes(pot.cdom); +pot.csize = sum(node_sizes(pot.cdom)); +pot.subtype = subtype; + +if nargin < 4 + can = cell(1, pot.dsize); + for i=1:pot.dsize + can{i} = cpot(cdom, node_sizes(cdom)); + end +end +pot.can = can; + +if nargin < 5 + mom = cell(1, pot.dsize); + for i=1:pot.dsize + mom{i} = mpot(cdom, node_sizes(cdom)); + end +end +pot.mom = mom; + +pot = class(pot, 'cgpot'); + diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/display.m new file mode 100644 index 00000000..78813bec --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/display.m @@ -0,0 +1,4 @@ +function display(pot) + +disp('conditional Gaussian potential object'); +disp(struct(pot)); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/divide_by_pot.m new file mode 100644 index 00000000..d3c6e7f9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/divide_by_pot.m @@ -0,0 +1,19 @@ +function bigpot = divide_by_pot(bigpot, smallpot) +% DIVIDE_BY_POT bigpot /= smallpot for cgpot +% bigpot = divide_by_pot(bigpot, smallpot) +% +% smallpot's domain must be a subset of bigpot's domain. + +bigpot = cg_mom_to_can(bigpot); +smallpot = cg_mom_to_can(smallpot); + +mask = find_equiv_posns(smallpot.ddom, bigpot.ddom); +for i=1:bigpot.dsize + if isempty(smallpot.ddom) + src = 1; + else + sub = ind2subv(bigpot.dsizes, i); + src = subv2ind(smallpot.dsizes, sub(mask)); + end + bigpot.can{i} = divide_by_pot(bigpot.can{i}, smallpot.can{src}); +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/domain_pot.m new file mode 100644 index 00000000..a09e39e6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/domain_pot.m @@ -0,0 +1,5 @@ +function dom = domain_pot(pot) +% DOMAIN_POT Return the domain of this cgpot. +% dom = domain_pot(pot) + +dom = pot.domain; diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_cts_evidence_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_cts_evidence_pot.m new file mode 100644 index 00000000..e5267bbd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_cts_evidence_pot.m @@ -0,0 +1,8 @@ +function pot = enter_cts_evidence_pot(pot, Y, y) +% function pot = enter_cts_evidence_pot(pot, Y, y) cgpot + + +pot = cg_mom_to_can(pot); +for i=1:pot.dsize + pot.can{i} = enter_cts_evidence_pot(pot.can{i}, Y, y); +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_discrete_evidence_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_discrete_evidence_pot.m new file mode 100644 index 00000000..972df966 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_discrete_evidence_pot.m @@ -0,0 +1,16 @@ +function pot = enter_discrete_evidence_pot(pot, Y, y) + +%ns = sparse(1, max(pot.domain)); +ns = zeros(1, max(pot.domain)); +ns(pot.ddom) = pot.dsizes; +ns(pot.cdom) = pot.csizes; + +ddom = pot.ddom; +S = prod(ns(ddom)); +sub = ind2subv(ns(ddom), 1:S); +mask = find_equiv_posns(Y, ddom); +sub(mask) = y; +ndx = subv2ind(ns(ddom), sub); + +pot.can = pot.can(ndx); +pot.mom = pot.mom(ndx); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/marginalize_pot.m new file mode 100644 index 00000000..4e666b1a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/marginalize_pot.m @@ -0,0 +1,153 @@ +function smallpot = marginalize_pot(bigpot, keep, maximize, useC) +% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, keep, maximize, useC) +% +% If maximize = 1, we raise an error. +% useC is ignored. + +if nargin < 3, maximize = 0; end +assert(~maximize); + + +sumover = mysetdiff(bigpot.domain, keep); +csumover = myintersect(sumover, bigpot.cdom); +dsumover = myintersect(sumover, bigpot.ddom); +dkeep = myintersect(keep, bigpot.ddom); +ckeep = myintersect(keep, bigpot.cdom); +%ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer +ns = zeros(1, max(bigpot.domain)); +ns(bigpot.ddom) = bigpot.dsizes; +ns(bigpot.cdom) = bigpot.csizes; + +% sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes. +% Similarly, prod(ns(dsumover))==1 is like isempty(dsumover) + +% Marginalize the cts parts. +% If we are in canonical form, we stay that way, since moment form might not exist. +% Besides, we would like to minimize the number of conversions. +if sum(ns(csumover)) > 0 + if bigpot.subtype == 'm' + for i=1:bigpot.dsize + bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep); + end + else + for i=1:bigpot.dsize + bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep); + end + end +end + +% If we are not marginalizing over any discrete nodes, we are done. +if prod(ns(dsumover))==1 + smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype); + return; +end + +% To marginalize the discrete parts, we partition the cts parts into those that depend +% on dkeep (i) and those that depend on on dsumover (j). + +I = prod(ns(dkeep)); +J = prod(ns(dsumover)); +C = sum(ns(ckeep)); +sum_map = find_equiv_posns(dsumover, bigpot.ddom); +keep_map = find_equiv_posns(dkeep, bigpot.ddom); +iv = zeros(1, length(bigpot.ddom)); % index vector + +% If in canonical form, marginalize if possible, else convert to moment form. +if 0 && bigpot.subtype == 'c' + p1 = zeros(I,J); + h1 = zeros(C,J,I); + K1 = zeros(C,C,J,I); + for i=1:I + keep_iv = ind2subv(ns(dkeep), i); + iv(keep_map) = keep_iv; + for j=1:J + sum_iv = ind2subv(ns(dsumover), j); + iv(sum_map) = sum_iv; + k = subv2ind(ns(bigpot.ddom), iv); + can = struct(bigpot.can{k}); % violate object privacy + p1(i,j) = exp(can.g); + if C > 0 % so mu1 and Sigma1 are non-empty + h1(:,j,i) = can.h; + K1(:,:,j,i) = can.K; + end + end + end + + % If the cts parts do not depend on j, we can just marginalize the weighting coefficient g. + jdepends = 0; + for i=1:I + for j=2:J + if ~approxeq(h1(:,j,i), h1(:,1,i)) | ~approxeq(K1(:,:,j,i), K1(:,:,1,i)) + jdepends = 1; + break + end + end + end + + if ~jdepends + %g2 = log(sum(p1, 2)); + g2 = zeros(I,1); + for i=1:I + s = sum(p1(i,:)); + if s > 0 + g2(i) = log(s); + end + end + h2 = h1; + K2 = K1; + can = cell(1,I); + j = 1; % arbitrary + for i=1:I + can{i} = cpot(ckeep, ns(ckeep), g2(i), h2(:,j,i), K2(:,:,j,i)); + end + smallpot = cgpot(dkeep, ckeep, ns, can, [], 'c'); + return; + else + % Since the cts parts depend on j, we must convert to moment form + bigpot = cg_can_to_mom(bigpot); + end +end + + +% Marginalize in moment form +bigpot = cg_can_to_mom(bigpot); + +% Now partition the moment components. +T1 = zeros(I,J); +mu1 = zeros(C,J,I); +Sigma1 = zeros(C,C,J,I); +for i=1:I + keep_iv = ind2subv(ns(dkeep), i); + iv(keep_map) = keep_iv; + for j=1:J + sum_iv = ind2subv(ns(dsumover), j); + iv(sum_map) = sum_iv; + k = subv2ind(ns(bigpot.ddom), iv); + mom = struct(bigpot.mom{k}); % violate object privacy + T1(i,j) = exp(mom.logp); + if C > 0 % so mu1 and Sigma1 are non-empty + mu1(:,j,i) = mom.mu; + Sigma1(:,:,j,i) = mom.Sigma; + end + end +end + +% Collapse the mixture of Gaussians +coef = mk_stochastic(T1); % coef must be convex combination +T2 = sum(T1,2); +T2 = T2 + (T2==0)*eps; +%if C > 0, disp('collapsing onto '); disp(leep); end +mu = []; +Sigma = []; +mom = cell(1,I); +for i=1:I + if C > 0 + [mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:)); + end + logp = log(T2(i)); + mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma); +end + +smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm'); + diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_by_pot.m new file mode 100644 index 00000000..0dcad3fc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_by_pot.m @@ -0,0 +1,19 @@ +function bigpot = multiply_by_pot(bigpot, smallpot, varargin) +% MULTIPLY_BY_POT bigpot *= smallpot for cgpot +% bigpot = multiply_by_pot(bigpot, smallpot) +% +% smallpot's domain must be a subset of bigpot's domain. + +bigpot = cg_mom_to_can(bigpot); +smallpot = cg_mom_to_can(smallpot); + +mask = find_equiv_posns(smallpot.ddom, bigpot.ddom); +for i=1:bigpot.dsize + if isempty(smallpot.ddom) + src = 1; + else + sub = ind2subv(bigpot.dsizes, i); + src = subv2ind(smallpot.dsizes, sub(mask)); + end + bigpot.can{i} = multiply_by_pot(bigpot.can{i}, smallpot.can{src}); +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_pots.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_pots.m new file mode 100644 index 00000000..8035d6de --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_pots.m @@ -0,0 +1,16 @@ +function T = multiply_pots(T1, T2) +% MULTIPLY_POTS Multiply a pair of dpots together pointwise (cgpot) +% T = multiply_pots(pots) + +ddom = myunion(T1.ddom, T2.ddom); +cdom = myunion(T1.cdom, T2.cdom); +dom = myunion(ddom, cdom); +ns = zeros(1, max(dom)); +ns(T1.ddom) = T1.dsizes; +ns(T2.ddom) = T2.dsizes; +ns(T1.cdom) = T1.csizes; +ns(T2.cdom) = T2.csizes; + +T = cgpot(ddom, cdom, ns); +T = multiply_by_pot(T, T1); +T = multiply_by_pot(T, T2); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/normalize_pot.m new file mode 100644 index 00000000..0c05ae5d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/normalize_pot.m @@ -0,0 +1,24 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the CG potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +% Marginalize down to [], so that the normalizing constant becomes Pr(E) +temp = cg_can_to_mom(marginalize_pot(pot, [])); +%loglik = temp.mom{1}.logp; +[temp2, loglik] = normalize_pot(temp.mom{1}); + +% Adjust scale factor to reflect the fact that the pot now represents Pr(X | E) instead of Pr(X,E). + +scale = -loglik; +if 1 +switch pot.subtype + case 'm' + for i=1:pot.dsize + pot.mom{i} = rescale_pot(pot.mom{i}, scale); + end + case 'c' + for i=1:pot.dsize + pot.can{i} = rescale_pot(pot.can{i}, scale); + end +end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/pot_to_marginal.m new file mode 100644 index 00000000..a651e192 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/pot_to_marginal.m @@ -0,0 +1,27 @@ +function m = pot_to_marginal(pot) +% POT_TO_MARGINAL Convert a cgpot to a marginal structure. +% m = pot_to_marginal(pot) + +pot = cg_can_to_mom(pot); +m.domain = pot.domain; +n = pot.csize; +d = length(pot.mom); +if n==0 + m.mu = []; + m.Sigma = []; +else + m.mu = zeros(n, d); + m.Sigma = zeros(n, n, d); +end +m.T = 0*myones(pot.dsizes); +for i=1:pot.dsize + s = struct(pot.mom{i}); % violate privacy of object + if n > 0 + m.mu(:,i) = s.mu; + m.Sigma(:,:,i) = s.Sigma; + end + m.T(i) = exp(s.logp); +end +if isvectorBNT(m.T) + m.T = m.T(:)'; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/set_domain_pot.m new file mode 100644 index 00000000..07364648 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/set_domain_pot.m @@ -0,0 +1,21 @@ +function pot = set_domain_pot(pot, domain) +% SET_DOMAIN_POT Change the domain of a potential (cgpot) +% pot = set_domain_pot(pot, domain) + +delta = domain(1) - pot.domain(1); +assert(all(domain == pot.domain + delta)); +pot.domain = pot.domain + delta; +pot.ddom = pot.ddom + delta; +pot.cdom = pot.cdom + delta; +cdomain = pot.cdom; +n = prod(pot.dsizes); +if(pot.subtype == 'm') + for i = 1: n + pot.mom{i} = set_domain_pot(pot.mom{i}, cdomain); + end +end +if(pot.subtype == 'c') + for i = 1: n + pot.can{i} = set_domain_pot(pot.can{i}, cdomain); + end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries new file mode 100644 index 00000000..7f4397ef --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries @@ -0,0 +1,14 @@ +/cpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/cpot_to_mpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/display.m/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/enter_cts_evidence_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/multiply_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/multiply_pots.m/1.1.1.1/Wed May 29 15:59:58 2002// +/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002// +/rescale_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/set_domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries.Log new file mode 100644 index 00000000..24f16336 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries.Log @@ -0,0 +1 @@ +A D/Old//// diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Repository new file mode 100644 index 00000000..12ecad25 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@cpot diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Entries new file mode 100644 index 00000000..3345f808 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Entries @@ -0,0 +1,3 @@ +/cpot_to_mpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/normalize_pot.convert.m/1.1.1.1/Wed May 29 15:59:58 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Repository new file mode 100644 index 00000000..bdb2b68d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@cpot/Old diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/cpot_to_mpot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/cpot_to_mpot.m new file mode 100644 index 00000000..05899b99 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/cpot_to_mpot.m @@ -0,0 +1,25 @@ +function mom = cpot_to_mpot(can) +% CPOT_TO_MPOT Convert a canonical potential to moment form. +% mom = cpot_to_mpot(can) + +[logp, mu, Sigma] = canonical_to_moment(can.g, can.h, can.K); +mom = mpot(can.domain, can.sizes, logp, mu, Sigma); + +%%%%%%% + +function [logp, mu, Sigma] = canonical_to_moment(g, h, K) +% CANONICAL_TO_MOMENT Convert canonical characteristics to moment form. +% [logp, mu, Sigma] = canonical_to_moment(g, h, K) + +if det(K)==0 + Sigma = inf*size(K); +else + Sigma = inv(K); +end +mu = Sigma*h; +n = length(mu); +if isempty(mu) + logp = g - 0.5*(log(det(K)) - n*log(2*pi)); +else + logp = g - 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu); +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/normalize_pot.convert.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/normalize_pot.convert.m new file mode 100644 index 00000000..59c473cc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/normalize_pot.convert.m @@ -0,0 +1,9 @@ +function [mom2, loglik] = normalize_pot(can) +% NORMALIZE_POT Convert the canonical potential Pr(X,E) into moment potential Pr(X|E) and return log Pr(E). +% [mom, loglik] = normalize_pot(can) + +mom = cpot_to_mpot(can); +mom = struct(mom); % violate privacy of object +loglik = mom.logp; +%mom.logp = 0; % now represents Pr(X | E) instead of Pr(X, E). +mom2 = mpot(mom.domain, mom.sizes, 0, mom.mu, mom.Sigma); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot.m new file mode 100644 index 00000000..b8dfbf03 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot.m @@ -0,0 +1,17 @@ +function pot = cpot(members, sizes, g, h, K) +% CPOT Make a canonical Gaussian potential. +% pot = cpot(members, sizes, g, h, K) +% +% All params default to 0 if omitted. + +n = sum(sizes); +if nargin < 3, g = 0; end +if nargin < 4, h = zeros(n,1); end +if nargin < 5, K = zeros(n,n); end + +pot.domain = members; +pot.sizes = sizes(:)'; +pot.g = g; +pot.h = h; +pot.K = K; +pot = class(pot, 'cpot'); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot_to_mpot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot_to_mpot.m new file mode 100644 index 00000000..73739d7e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot_to_mpot.m @@ -0,0 +1,29 @@ +function mom = cpot_to_mpot(can) +% CPOT_TO_MPOT Convert a canonical potential to moment form. +% mom = cpot_to_mpot(can) + +[logp, mu, Sigma] = canonical_to_moment(can.g, can.h, can.K); +mom = mpot(can.domain, can.sizes, logp, mu, Sigma); + +%%%%%%% + +function [logp, mu, Sigma] = canonical_to_moment(g, h, K) +% CANONICAL_TO_MOMENT Convert canonical characteristics to moment form. +% [logp, mu, Sigma] = canonical_to_moment(g, h, K) + +n = length(K); +if isempty(K) + logp = g - 0.5*(log(det(K)) - n*log(2*pi)); + Sigma = []; + mu = []; +else + if det(K)==0 + Sigma = inf*ones(n,n); + mu = zeros(n,1); % if the precision is zero, the mean is arbitrary + logp = g; % the scaling factor for the uniform distribution is 1 + else + Sigma = inv(K); + mu = Sigma*h; + logp = g - 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu); + end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/display.m new file mode 100644 index 00000000..35a4f0b3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/display.m @@ -0,0 +1,4 @@ +function display(pot) + +disp('canonical potential object'); +disp(struct(pot)); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/divide_by_pot.m new file mode 100644 index 00000000..1e2bb2b9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/divide_by_pot.m @@ -0,0 +1,13 @@ +function bigpot = divide_by_pot(bigpot, smallpot) +% DIVIDE_BY_POT bigpot /= smallpot for cpot +% bigpot = divide_by_pot(bigpot, smallpot) +% +% smallpot's domain must be a subset of bigpot's domain. + +bigpot.g = bigpot.g - smallpot.g; +if sum(smallpot.sizes) > 0 + mask = find_equiv_posns(smallpot.domain, bigpot.domain); + u = block(mask, bigpot.sizes); + bigpot.h(u) = bigpot.h(u) - smallpot.h; + bigpot.K(u, u) = bigpot.K(u, u) - smallpot.K; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/domain_pot.m new file mode 100644 index 00000000..aa5ce392 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/domain_pot.m @@ -0,0 +1,5 @@ +function dom = domain_pot(pot) +% DOMAIN_POT Return the domain of this cpot. +% dom = domain_pot(pot) + +dom = pot.domain; diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/enter_cts_evidence_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/enter_cts_evidence_pot.m new file mode 100644 index 00000000..246729ae --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/enter_cts_evidence_pot.m @@ -0,0 +1,15 @@ +function pot = enter_cts_evidence_pot(pot, Y, y) +% function pot = enter_cts_evidence_pot(pot, Y, y) (cpot) + +ns = sparse(1, max(pot.domain)); +ns(pot.domain) = pot.sizes; + +X = mysetdiff(pot.domain, Y); +[hx, hy, KXX, KXY, KYX, KYY] = partition_matrix_vec(pot.h, pot.K, X, Y, ns); +pot.g = pot.g + hy'*y - 0.5*y'*KYY*y; +if ~isempty(X) + pot.h = hx - KXY*y; + pot.K = KXX; +end + +pot.sizes(find_equiv_posns(Y,pot.domain)) = 0; diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/marginalize_pot.m new file mode 100644 index 00000000..ae9ee839 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/marginalize_pot.m @@ -0,0 +1,31 @@ +function smallpot = marginalize_pot(bigpot, keep, maximize, useC) +% MARGINALIZE_POT Marginalize a cpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, keep, maximize, useC) +% +% The maximize argument is ignored - maxing out a Gaussian is the same as summing it out, +% since the mode and mean are equal. +% The useC argument is ignored. + +node_sizes = sparse(1, max(bigpot.domain)); +node_sizes(bigpot.domain) = bigpot.sizes; +sum_over = mysetdiff(bigpot.domain, keep); + +if sum(node_sizes(sum_over))==0 % isempty(sum_over) + %smallpot = bigpot; + smallpot = cpot(keep, node_sizes(keep), bigpot.g, bigpot.h, bigpot.K); +else + [h1, h2, K11, K12, K21, K22] = partition_matrix_vec(bigpot.h, bigpot.K, sum_over, keep, node_sizes); + n = length(h1); + K11inv = inv(K11); + g = bigpot.g + 0.5*(n*log(2*pi) - log(det(K11)) + h1'*K11inv*h1); + if length(h2) > 0 % ~isempty(keep) % we are are actually keeping something + A = K21*K11inv; + h = h2 - A*h1; + K = K22 - A*K12; + else + h = []; + K = []; + end + smallpot = cpot(keep, node_sizes(keep), g, h, K); +end + diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_by_pot.m new file mode 100644 index 00000000..1ab24f55 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_by_pot.m @@ -0,0 +1,13 @@ +function bigpot = multiply_by_pot(bigpot, smallpot, varargin) +% MULTIPLY_BY_POT bigpot *= smallpot for cpot +% bigpot = multiply_by_pot(bigpot, smallpot) +% +% smallpot's domain must be a subset of bigpot's domain. + +bigpot.g = bigpot.g + smallpot.g; +if sum(smallpot.sizes) > 0 + mask = find_equiv_posns(smallpot.domain, bigpot.domain); + u = block(mask, bigpot.sizes); + bigpot.h(u) = bigpot.h(u) + smallpot.h; + bigpot.K(u, u) = bigpot.K(u, u) + smallpot.K; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_pots.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_pots.m new file mode 100644 index 00000000..3506e36e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_pots.m @@ -0,0 +1,11 @@ +function T = multiply_pots(T1, T2) +% MULTIPLY_POTS Multiply a pair of dpots together pointwise (cpot) +% T = multiply_pots(pots) + +dom = myunion(T1.domain, T2.domain); +ns = sparse(1, max(dom)); +ns(T1.domain) = T1.sizes; +ns(T2.domain) = T2.sizes; +T = cpot(dom, ns(dom)); +T = multiply_by_pot(T, T1); +T = multiply_by_pot(T, T2); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/normalize_pot.m new file mode 100644 index 00000000..c6f89dc3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/normalize_pot.m @@ -0,0 +1,10 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the canonical potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +mom = cpot_to_mpot(pot); % move the normalizing constant out of g, to reveal the coefficient +%loglik = scaling_factor_pot(mom); +%loglik = mom.logp; +[temp, loglik] = normalize_pot(mom); +pot.g = pot.g - loglik; + diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/pot_to_marginal.m new file mode 100644 index 00000000..4c3fb79f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/pot_to_marginal.m @@ -0,0 +1,6 @@ +function m = pot_to_marginal(pot) +% POT_TO_MARGINAL Convert a cpot to a marginal structure. +% m = pot_to_marginal(pot) + +mom = cpot_to_mpot(pot); +m = pot_to_marginal(mom); diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/rescale_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/rescale_pot.m new file mode 100644 index 00000000..ad962455 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/rescale_pot.m @@ -0,0 +1,5 @@ +function pot = rescale_pot(pot, s) +% RESCALE_POT Add a constant to the cpot scale factor. +% pot = rescale_pot(pot, s) + +pot.g = pot.g + s; diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/set_domain_pot.m new file mode 100644 index 00000000..2178c51d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/set_domain_pot.m @@ -0,0 +1,5 @@ +function pot = set_domain_pot(pot, domain) +% SET_DOMAIN_POT Change the domain of a potential (dpot) +% pot = set_domain_pot(pot, domain) + +pot.domain = domain; diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Entries new file mode 100644 index 00000000..d37ad7eb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Entries @@ -0,0 +1,17 @@ +/approxeq_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/display.m/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_pot.m/1.1.1.1/Thu Aug 5 15:25:08 2004// +/domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/dpot.m/1.1.1.1/Tue Oct 1 19:04:44 2002// +/dpot_to_table.m/1.1.1.1/Wed May 29 15:59:58 2002// +/find_most_prob_entry.m/1.1.1.1/Sun Jun 16 19:06:20 2002// +/get_fields.m/1.1.1.1/Wed May 29 15:59:58 2002// +/marginalize_pot.m/1.1.1.1/Wed Aug 4 19:59:14 2004// +/multiply_by_pot.m/1.1.1.1/Wed Aug 4 19:59:14 2004// +/multiply_pots.m/1.1.1.1/Wed May 29 15:59:58 2002// +/normalize_pot.m/1.1.1.1/Wed Aug 4 15:54:48 2004// +/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002// +/set_domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/subsasgn.m/1.1.1.1/Wed Apr 27 18:34:48 2005// +/subsref.m/1.1.1.1/Wed Apr 27 18:34:48 2005// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Repository new file mode 100644 index 00000000..485f27df --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@dpot diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/approxeq_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/approxeq_pot.m new file mode 100644 index 00000000..b0f8d0ed --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/approxeq_pot.m @@ -0,0 +1,5 @@ +function p = approxeq_pot(A, B, tol) + +if nargin < 3, tol = 1e-3; end + +p = approxeq(A.T, B.T, tol); diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/display.m new file mode 100644 index 00000000..1b0647bb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/display.m @@ -0,0 +1,4 @@ +function display(pot) + +disp('discrete potential object'); +disp(struct(pot)); diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/divide_by_pot.m new file mode 100644 index 00000000..aae4790c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/divide_by_pot.m @@ -0,0 +1,14 @@ +function Tbig = divide_by_pot(Tbig, Tsmall) +% DIVIDE_BY_POT Tbig /= Tsmall +% Tbig = divide_by_pot(Tbig, Tsmall) +% +% Tsmall's domain must be a subset of Tbig's domain. + +%process sparse dpot, we do not concern only one of the two pots is sparse +if issparse(Tbig.T) && issparse(Tsmall.T) + Tbig.T = divide_by_sparse_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes); +else + Tbig.T = divide_by_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes); +end + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/domain_pot.m new file mode 100644 index 00000000..a67cef02 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/domain_pot.m @@ -0,0 +1,5 @@ +function dom = domain_pot(pot) +% DOMAIN_POT Return the domain of this dpot. +% dom = domain_pot(pot) + +dom = pot.domain; diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot.m new file mode 100644 index 00000000..74b5dd9e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot.m @@ -0,0 +1,26 @@ +function pot = dpot(domain, sizes, T) +% DPOT Make a discrete (sparse) potential. +% pot = dpot(domain, sizes, T, spar) +% +% sizes(i) is the size of the i'th domain element. +% T defaults to all 1s. + +%assert(length(sizes) == length(domain)); + +pot.domain = domain(:)'; % so we can see it when we display +if nargin < 3 + pot.T = myones(sizes); + %pot.T = ones(1,prod(sizes)); % 1D vector +else + if isempty(T) + pot.T = []; + else + if issparse(T) + pot.T = T; + else + pot.T = myreshape(T, sizes); + end + end +end +pot.sizes = sizes(:)'; +pot = class(pot, 'dpot'); diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot_to_table.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot_to_table.m new file mode 100644 index 00000000..bd2d3d98 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot_to_table.m @@ -0,0 +1,3 @@ +function T = dpot_to_table(pot) + +T = pot.T; diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/find_most_prob_entry.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/find_most_prob_entry.m new file mode 100644 index 00000000..d30561a8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/find_most_prob_entry.m @@ -0,0 +1,10 @@ +function [indices, pot] = find_most_prob_entry(pot) +% function [indices, pot] = find_most_prob_entry(pot) +% function [indices, pot] = find_most_prob_entry(pot) +% Find the indices of the argmax, and set all other enties to 0. + +%indices = argmax(pot.T); +[m i] = max(pot.T(:)); +indices = ind2subv(pot.sizes, i); +pot.T = 0*myones(pot.sizes); +pot.T(i) = m; diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/get_fields.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/get_fields.m new file mode 100644 index 00000000..2400bc99 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/get_fields.m @@ -0,0 +1,13 @@ +function val = get_params(pot, name) +% GET_PARAMS Accessor function for a field (dpot) +% val = get_params(pot, name) +% +% e.g., get_params(pot, 'table') or 'domain' + +switch name + case 'table', val = pot.T; + case 'domain', val = pot.domain; + otherwise, + error(['invalid field name ' name]); +end + diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/marginalize_pot.m new file mode 100644 index 00000000..d38c5c82 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/marginalize_pot.m @@ -0,0 +1,17 @@ +function smallpot = marginalize_pot(bigpot, onto, maximize) +% MARGINALIZE_POT Marginalize a dpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, onto, maximize) +% +% 'onto' must be in ascending order. + +if nargin < 3, maximize = 0; end + +ns = zeros(1, max(bigpot.domain)); +ns(bigpot.domain) = bigpot.sizes; +%assert(isequal(bigpot.sizes, mysize(bigpot.T))); % may fail if there are trailing dimensions of size 1 +if issparse(bigpot.T) + smallT = marg_sparse_table(bigpot.T, bigpot.domain, bigpot.sizes, onto, maximize); +else + smallT = marg_table(bigpot.T, bigpot.domain, bigpot.sizes, onto, maximize); +end +smallpot = dpot(onto, ns(onto), smallT); diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_by_pot.m new file mode 100644 index 00000000..ad7caaa7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_by_pot.m @@ -0,0 +1,13 @@ +function Tbig = multiply_by_pot(Tbig, Tsmall) +% MULTIPLY_BY_POT Tbig *= Tsmall +% Tbig = multiply_by_pot(Tbig, Tsmall) +% +% Tsmall's domain must be a subset of Tbig's domain. + +%process sparse dpot, we do not consider only one of the two pots is sparse +if issparse(Tbig.T) && issparse(Tsmall.T) + Tbig.T = mult_by_sparse_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes); +else + Tbig.T = mult_by_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes); +end + diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_pots.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_pots.m new file mode 100644 index 00000000..0188f01b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_pots.m @@ -0,0 +1,12 @@ +function T = multiply_pots(T1, T2) +% MULTIPLY_POTS Multiply a pair of dpots together pointwise. +% T = multiply_pots(pots) + +dom = myunion(T1.domain, T2.domain); +%ns = sparse(1, max(dom)); % causes problems in myreshape on NT +ns = zeros(1, max(dom)); +ns(T1.domain) = T1.sizes; +ns(T2.domain) = T2.sizes; +T = dpot(dom, ns(dom)); +T = multiply_by_pot(T, T1); +T = multiply_by_pot(T, T2); diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/normalize_pot.m new file mode 100644 index 00000000..1bc43cdb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/normalize_pot.m @@ -0,0 +1,12 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the discrete potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +if isempty(pot.T) %add to process sparse + loglik = 0; + return; +end +[pot.T, lik] = normalise(pot.T); +loglik = log(lik + (lik==0)*eps); + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/pot_to_marginal.m new file mode 100644 index 00000000..36ed70ea --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/pot_to_marginal.m @@ -0,0 +1,12 @@ +function m = pot_to_marginal(pot) +% POT_TO_MARGINAL Convert a dpot to a marginal structure. +% m = pot_to_marginal(pot) + +m.domain = pot.domain; +m.T = pot.T; +m.mu = []; +m.Sigma = []; + +%if isvector(m.T) +% m.T = m.T(:); +%end diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/set_domain_pot.m new file mode 100644 index 00000000..2178c51d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/set_domain_pot.m @@ -0,0 +1,5 @@ +function pot = set_domain_pot(pot, domain) +% SET_DOMAIN_POT Change the domain of a potential (dpot) +% pot = set_domain_pot(pot, domain) + +pot.domain = domain; diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/subsasgn.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsasgn.m new file mode 100644 index 00000000..90d43d70 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsasgn.m @@ -0,0 +1,3 @@ +function B = subsasgn(A, S, B) + +B = builtin('subsasgn', A, S, B); diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/subsref.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsref.m new file mode 100644 index 00000000..4a48bffe --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsref.m @@ -0,0 +1,3 @@ +function B = subsref(A, S) + +B = builtin('subsref', A, S); diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Entries new file mode 100644 index 00000000..c542af05 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Entries @@ -0,0 +1,9 @@ +/display.m/1.1.1.1/Wed May 29 15:59:58 2002// +/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/mpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/mpot_to_cpot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002// +/rescale_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/set_domain_pot.m/1.1.1.1/Wed Jul 30 13:37:52 2003// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Repository new file mode 100644 index 00000000..1d9311e4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@mpot diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/display.m new file mode 100644 index 00000000..90f0e2fc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/display.m @@ -0,0 +1,4 @@ +function display(pot) + +disp('moment Gaussian potential object'); +disp(struct(pot)); diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/marginalize_pot.m new file mode 100644 index 00000000..88751c62 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/marginalize_pot.m @@ -0,0 +1,27 @@ +function smallpot = marginalize_pot(bigpot, keep, maximize, useC) +% MARGINALIZE_POT Marginalize a mpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, keep, maximize, useC) +% +% The maximize argument is ignored - maxing out a Gaussian is the same as summing it out, +% since the mode and mean are equal. +% The useC argument is ignored. + + +node_sizes = sparse(1, max(bigpot.domain)); +node_sizes(bigpot.domain) = bigpot.sizes; +sum_over = mysetdiff(bigpot.domain, keep); + +[logp, mu, Sigma] = marginalize_gaussian(bigpot.logp, bigpot.mu, bigpot.Sigma, ... + keep, sum_over, node_sizes); +smallpot = mpot(keep, node_sizes(keep), logp, mu, Sigma); + +%%%%%% + +function [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns) +% MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian. +% [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns) +% +% sizes(i) is the size of the i'th block in domain. + +[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns); +logpX = logp; % Lauritzen (1996) p161 diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot.m new file mode 100644 index 00000000..1c790f8a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot.m @@ -0,0 +1,17 @@ +function pot = mpot(members, sizes, logp, mu, Sigma) +% MPOT Make a moment Gaussian potential. +% pot = mpot(members, sizes, logp, mu, Sigma) +% +% All params default to 0 if omitted. + +n = sum(sizes); +if nargin < 3, logp = 0; end +if nargin < 4, mu = zeros(n,1); end +if nargin < 5, Sigma = zeros(n,n); end + +pot.domain = members; +pot.sizes = sizes; +pot.logp = logp; +pot.mu = mu; +pot.Sigma = Sigma;zeros(n,n); +pot = class(pot, 'mpot'); diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot_to_cpot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot_to_cpot.m new file mode 100644 index 00000000..ffb3192b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot_to_cpot.m @@ -0,0 +1,21 @@ +function can = mpot_to_cpot(mom) +% MPOT_TO_CPOT Convert a moment potential to canonical form. +% mom = mpot_to_cpot(can) + +[g, h, K] = moment_to_canonical(mom.logp, mom.mu, mom.Sigma); +can = cpot(mom.domain, mom.sizes, g, h, K); + +%%%%%%%%%%% + +function [g, h, K] = moment_to_canonical(logp, mu, Sigma) +% MOMENT_TO_CANONICAL Convert moment characteristics to canonical form. +% [g, h, K] = moment_to_canonical(logp, mu, Sigma) + +K = inv(Sigma); +h = K*mu; +n = length(K); +if isempty(mu) + g = logp + 0.5*(log(det(K)) - n*log(2*pi)); +else + g = logp + 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu); +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/normalize_pot.m new file mode 100644 index 00000000..a9c2e935 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/normalize_pot.m @@ -0,0 +1,6 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the moment potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +loglik = pot.logp; +pot.logp = 0; diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/pot_to_marginal.m new file mode 100644 index 00000000..b5bd041e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/pot_to_marginal.m @@ -0,0 +1,13 @@ +function m = pot_to_marginal(pot) +% POT_TO_MARGINAL Convert a mpot to a marginal structure. +% m = pot_to_marginal(pot) + +m.domain = pot.domain; +m.T = exp(pot.logp); +m.mu = pot.mu; +m.Sigma = pot.Sigma; + +if isvectorBNT(m.T) + m.T = m.T(:)'; +end + diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/rescale_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/rescale_pot.m new file mode 100644 index 00000000..39d1809d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/rescale_pot.m @@ -0,0 +1,5 @@ +function pot = rescale_pot(pot, s) +% RESCALE_POT Add a constant to the mpot scale factor. +% pot = rescale_pot(pot, s) + +pot.logp = pot.logp + s; diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/set_domain_pot.m new file mode 100644 index 00000000..91d5a7bb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/set_domain_pot.m @@ -0,0 +1,5 @@ +function pot = set_domain_pot(pot, domain) +% SET_DOMAIN_POT Change the domain of a potential (mpot) +% pot = set_domain_pot(pot, domain) + +pot.domain = domain; diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Entries new file mode 100644 index 00000000..d5b4bd91 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Entries @@ -0,0 +1,6 @@ +/marginalize_pot.m/1.1.1.1/Sun May 19 22:11:08 2002// +/normalize_pot.m/1.1.1.1/Sun May 19 22:11:08 2002// +/reduce_pot.m/1.1.1.1/Tue Mar 11 17:37:02 2003// +/rescale_pot.m/1.1.1.1/Sun May 19 22:11:08 2002// +/scgcpot.m/1.1.1.1/Sun May 19 22:11:08 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Repository new file mode 100644 index 00000000..18dc4d80 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@scgcpot diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/marginalize_pot.m new file mode 100644 index 00000000..c4ca9b68 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/marginalize_pot.m @@ -0,0 +1,7 @@ +function smallpot = marginalize_pot(bigpot, keepdom, sumoverdom, nodesizes) +% MARGINALIZE_POT Marginalize a mpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, keep) + +keepsize = sum(nodesizes(keepdom)); +[A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(bigpot.A, bigpot.B, bigpot.C, keepdom, sumoverdom, nodesizes); +smallpot = scgcpot(keepsize, bigpot.ctailsize, bigpot.p, A1, B1, C11); \ No newline at end of file diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/normalize_pot.m new file mode 100644 index 00000000..8900c66e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/normalize_pot.m @@ -0,0 +1,6 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the element of stable conditional gaussian potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +loglik = log(pot.p); +pot.p = 1; diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/reduce_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/reduce_pot.m new file mode 100644 index 00000000..466a7a72 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/reduce_pot.m @@ -0,0 +1,32 @@ +function [reduced_pot,successful] = reduce(pot,tailnodes) +% Executes the reduce operation defined in +% Stable Local Computation with Conditional Gaussian Distributions +% Steffen L. Lauritzen +% Frank Jensen +% September 1999 +% The potential pot is reduced if B contains any zero columns +% The test are restricted to the positions in tailnodes. +% Any columns successfully deleted are entered in the array successful + +if nargin < 2 + tailnodes = 1:pot.ctailsize; +end + +successful = []; + +% Look for all columns beeing equal to zero +for i = tailnodes + if ~any(pot.B(:,i)) + successful = [successful i]; + end +end + +remain = mysetdiff(1:pot.ctailsize,successful); + +% Erase the zero-columns and decrease the tailsize +pot.B = pot.B(:,remain); +pot.ctailsize = pot.ctailsize - length(successful); + +% Return the reduced potential +reduced_pot = pot; + diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/rescale_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/rescale_pot.m new file mode 100644 index 00000000..996a6351 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/rescale_pot.m @@ -0,0 +1,5 @@ +function pot = rescale_pot(pot, s) +% RESCALE_POT Add a constant to the mpot scale factor. +% pot = rescale_pot(pot, s) + +pot.p = pot.p*s; diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/scgcpot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/scgcpot.m new file mode 100644 index 00000000..6d53a16f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/scgcpot.m @@ -0,0 +1,49 @@ +function pot = scgcpot(cheadsize, ctailsize, p, A, B, C) +% SCGCPOT Make a base object of stable conditional gaussian potential. +% pot = scgcpot(cheadsize, ctailsize, p, A, B, C) +% +% cheadsize is the demension of head nodes. +% ctailsize is the demension of tail nodes. +% r = cheadsize, s = ctailsize +% p is discrete probability. +% A is table of r*1 vectors; +% B is r*s matrices +% C is r*r positive semidefinite symmetric matrices + +if nargin < 3 + p = 1; +end +if nargin < 4 + A = zeros(cheadsize,1); +end +if nargin < 5 + B = zeros(cheadsize,ctailsize); +end +if nargin < 6 + C = zeros(cheadsize,cheadsize); +end + +if isempty(A) + A = zeros(cheadsize,1); +end +if isempty(B) + B = zeros(cheadsize,ctailsize); +end +if isempty(C) + C = zeros(cheadsize,cheadsize); +end + +pot.cheadsize = cheadsize; +pot.ctailsize = ctailsize; + +pot.p = p; +pot.A = A; +pot.B = B; +pot.C = C; +%if cheadsize == 0 +% pot.A = []; +%end +%if ctailsize == 0 +% pot.B = []; +%end +pot = class(pot, 'scgcpot'); diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Entries new file mode 100644 index 00000000..c990d57d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Entries @@ -0,0 +1,12 @@ +/README/1.1.1.1/Thu Mar 20 15:07:16 2003// +/combine_pots.m/1.1.1.1/Tue Mar 11 17:49:28 2003// +/complement_pot.m/1.1.1.1/Wed May 21 13:49:34 2003// +/direct_combine_pots.m/1.1.1.1/Sun May 19 22:11:08 2002// +/extension_pot.m/1.1.1.1/Fri Jan 24 12:52:34 2003// +/marginalize_pot.m/1.1.1.1/Tue Mar 11 17:06:08 2003// +/normalize_pot.m/1.1.1.1/Wed May 21 13:49:44 2003// +/pot_to_marginal.m/1.1.1.1/Sun May 19 22:11:08 2002// +/recursive_combine_pots.m/1.1.1.1/Wed May 21 13:49:48 2003// +/reduce_pot.m/1.1.1.1/Tue Mar 11 18:07:12 2003// +/scgpot.m/1.1.1.1/Tue Mar 11 14:04:48 2003// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Repository new file mode 100644 index 00000000..ffef4774 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@scgpot diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/README b/sourcecodes/bnt-master/BNT/potentials/@scgpot/README new file mode 100644 index 00000000..1a3c9d16 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/README @@ -0,0 +1,11 @@ +% Stable conditional Gaussian inference +% Written by Rainer Deventer + + +@techreport{Lauritzen99, + author = "S. Lauritzen and F. Jensen", + title = "Stable Local Computation with Conditional {G}aussian Distributions", + year = 1999, + number = "R-99-2014", + institution = "Dept. Math. Sciences, Aalborg Univ." +} diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/combine_pots.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/combine_pots.m new file mode 100644 index 00000000..2315176a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/combine_pots.m @@ -0,0 +1,24 @@ +function pot = combine_pots(pot1, pot2) +% COMBINE_POTS combine two potentials +% pot = combine_pots(pot1, pot2) + +% Reduce both potentials before trying to combine them. +% Cf. "Stable Local computation with Conditional Gaussian Distributions", page 9 +% Consider again two potentials with minimal tail + +% Guarantee minimal tails. If pot1 or pot2 are minimal, they are not changed +pot1 = reduce_pot(pot1); +pot2 = reduce_pot(pot2); + +%if the intersect set of these two potentials' head conts. combination is undifined +if ~isempty( myintersect(pot1.cheaddom, pot2.cheaddom) ) + return; +end + +if isempty( myintersect(pot1.domain, pot2.cheaddom) ) | isempty( myintersect(pot2.domain, pot1.cheaddom)) + % if satisfy the condition of directed combine + pot = direct_combine_pots(pot1, pot2); +else + % perform recursive combine + pot = recursive_combine_pots(pot1, pot2); +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/complement_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/complement_pot.m new file mode 100644 index 00000000..25f1c1c6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/complement_pot.m @@ -0,0 +1,221 @@ +function [margpot, comppot] = complement_pot(pot, keep) +% COMPLEMENT_POT complement means decompose of a potential into its strong marginal and +% its complement corresponds exactly to the decomposition of a probability distribution +% into its marginal and conditional +% [margpot, comppot] = complement_pot(pot, keep) + +% keep can only include continuous head nodes and discrete nodes +% margpot is the stable CG potential of keep nodes +% comppot is the stable CG potential of others in corresponds exactly to +% the discomposition of a probability distribution of its marginal and conditional + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Calculation of the marginal requires integration over % +% all variables in csumover. Thus cheadkeep contains all % +% continuous variables in the marginal potential % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%keyboard; +csumover = mysetdiff(pot.cheaddom, keep); +cheadkeep = mysetdiff(pot.cheaddom, csumover); + +nodesizes = zeros(1, max(pot.domain)); +nodesizes(pot.ddom) = pot.dsizes; +nodesizes(pot.cheaddom) = pot.cheadsizes; +nodesizes(pot.ctaildom) = pot.ctailsizes; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Description of the variables in the marginal domain % +% For the calculation of a strong marginal first integration % +% over all continuous variables in the head takes place. % +% The calculation of the marginal over the head variables % +% might result in a smaller or empty tail % +% If there are no head variables, and therefore no tail % +% variables, left marginalisation over discrete variables % +% may take place % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +margdom = mysetdiff(pot.domain,keep); +% margddom = pot.ddom; +margcheaddom = cheadkeep; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Marginalisation over discrete variables is only allowed when % +% the tail is empty % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +margddom = myintersect(pot.ddom,keep); % Discrete domain of marginal +margctaildom = myintersect(pot.ctaildom,keep); % Tail domain +assert(isempty(mysetdiff(pot.ddom,margddom)) | isempty(margctaildom)) + + +%margctaildom = pot.ctaildom; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Even if marginalisation over continuous variables is only defined % +% for head variables, the marginalisation over haed-variables might % +% result in a smaller tail % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +margctaildom = myintersect(pot.ctaildom,keep); + +margcheadsizes = nodesizes(margcheaddom); +margcheadsize = sum(margcheadsizes); +margctailsizes = nodesizes(margctaildom); +margctailsize = sum(margctailsizes); + +compdom = pot.domain; +compddom = pot.ddom; +compcheaddom = csumover; +compctaildom = myunion(pot.ctaildom, cheadkeep); +compcheadsizes = nodesizes(compcheaddom); +compcheadsize = sum(compcheadsizes); +compctailsizes = nodesizes(compctaildom); +compctailsize = sum(compctailsizes); + +dkeep = myintersect(pot.ddom, keep); +%if dom is only contain discrete node +if isempty(pot.cheaddom) + dsumover = mysetdiff(pot.ddom, dkeep); + + if isempty(dsumover) + margpot = pot; + comppot = scgpot([], [], [], []); + return; + end + + + I = prod(nodesizes(dkeep)); + J = prod(nodesizes(dsumover)); + sum_map = find_equiv_posns(dsumover, pot.ddom); + keep_map = find_equiv_posns(dkeep, pot.ddom); + iv = zeros(1, length(pot.ddom)); % index vector + p1 = zeros(I,J); + for i=1:I + keep_iv = ind2subv(nodesizes(dkeep), i); + iv(keep_map) = keep_iv; + for j=1:J + sum_iv = ind2subv(nodesizes(dsumover), j); + iv(sum_map) = sum_iv; + k = subv2ind(nodesizes(pot.ddom), iv); + potc = struct(pot.scgpotc{k}); % violate object privacy + p1(i,j) = potc.p; + end + end + p2 = sum(p1,2); + p2 = p2 + (p2==0)*eps; + + margscpot = cell(1, I); + compscpot = cell(1, I*J); + iv = zeros(1, length(pot.ddom)); % index vector + for i=1:I + margscpot{i} = scgcpot(0, 0, p2(i)); + keep_iv = ind2subv(nodesizes(dkeep), i); + iv(keep_map) = keep_iv; + for j=1:J + sum_iv = ind2subv(nodesizes(dsumover), j); + iv(sum_map) = sum_iv; + k = subv2ind(nodesizes(pot.ddom), iv); + q = p1(i,j)/p2(i); + compscpot{k} = scgcpot(0, 0, q); + end + end + + margpot = scgpot(dkeep, [], [], nodesizes, margscpot); + comppot = scgpot(pot.ddom, [], [], nodesizes,compscpot); + return; +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% head of the potential is not empty % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +dsize = pot.dsize; +compscpot = cell(1, dsize); + +fmaskh = find_equiv_posns(margcheaddom, compctaildom); +fmaskt = find_equiv_posns(margctaildom, compctaildom); + +fh = block(fmaskh, compctailsizes); +ft = block(fmaskt, compctailsizes); + + +if ~isempty(margcheaddom) + for i=1:dsize + potc = struct(pot.scgpotc{i}); + q = 1; + p = potc.p; + [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(potc.A, potc.B, potc.C, margcheaddom, compcheaddom, nodesizes); + + if ~isempty(margcheaddom) + margscpot{i} = scgcpot(margcheadsize, margctailsize, p, A1, B1, C11); + else + margscpot{i} = scgcpot(margcheadsize, margctailsize, p); + end + + if ~isempty(compcheaddom) + if ~isempty(margcheaddom) + E = A2 - C21*pinv(C11)*A1; + tmp1 = C21*pinv(C11); + tmp2 = B2 - C21*pinv(C11)*B1; + F = zeros(compcheadsize, compctailsize); + F(:, fh) = tmp1; + F(:, ft) = tmp2; + G = C22 - C21*pinv(C11)*C12; + else + E = A2; + F = B2; + G = C22; + end + compscpot{i} = scgcpot(compcheadsize, compctailsize, q, E, F, G); + else + compscpot{i} = scgcpot(compcheadsize, 0, q); + end + if isempty(margcheaddom) + margpot = scgpot(margddom, [], [], nodesizes, margscpot); + else + margpot = scgpot(margddom, margcheaddom, margctaildom, nodesizes, margscpot); + end + end +else + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Marginalisation took place over all head variables. % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Calculate the strong marginal % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + margpot = marginalize_pot(pot,keep); + mPot = struct(margpot); + for i =1:dsize + potc = struct(pot.scgpotc{i}); + % Get the probability of the original potential % + q = potc.p; + + % Get the configuration defined by the index i% + config = ind2subv(pot.dsizes,i); + + % Calculate the corresponding configuration in the marginal potential + if isempty(margpot.dsizes) + % keep == [] + indMargPot = 1; + else + equivPos = find_equiv_posns(dkeep,pot.ddom); + indMargPot = subv2ind(margpot.dsizes,config(equivPos)); + end + % Figure out the corresponding marginal potential + mPotC = struct(mPot.scgpotc{indMargPot}); + p = mPotC.p; + if p == 0 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % The following assignment is correct as p is only zero if q is also zero % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + compscpot{i} = scgcpot(compcheadsize,compctailsize,0,potc.A,potc.B,potc.C); + else + compscpot{i} = scgcpot(compcheadsize,compctailsize,q/p,potc.A,potc.B,potc.C); + end + end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Put all components in one potential % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if isempty(compcheaddom) + comppot = scgpot(compddom, [], [], nodesizes,compscpot); +else + comppot = scgpot(compddom, compcheaddom, compctaildom, nodesizes,compscpot); +end + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/direct_combine_pots.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/direct_combine_pots.m new file mode 100644 index 00000000..deeb002d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/direct_combine_pots.m @@ -0,0 +1,161 @@ +function pot = direct_combine_pots(pot1, pot2) +% DIRECTED_COMBINE_POTS The combination operation corresponds to ordinary composition of conditional distributions. +% In some sense is similar to that of forming disjoint union of set. +% pot = direct_combine_pots(pot1, pot2) + +% directed combine can be performed under the conditon that the head node set of pot1 is disjoint from the domain of +% pot2 or vice versa. if the last conditon was satisfied we exchange the pot1 and pot2 firstly then perform the operation. +% If neither of them was satified the directed combine is undifined. + + +if isempty( myintersect(pot1.domain, pot2.cheaddom) ) + pot1 = pot1; + pot2 = pot2; +elseif isempty( myintersect(pot2.domain, pot1.cheaddom)) + temppot = pot1; + pot1 = pot2; + pot2 = temppot; +else + assert(0); + return; +end + +domain = myunion(pot1.domain, pot2.domain); +nodesizes = zeros(1,max(domain)); +nodesizes(pot2.ctaildom) = pot2.ctailsizes; +nodesizes(pot2.cheaddom) = pot2.cheadsizes; +nodesizes(pot2.ddom) = pot2.dsizes; +nodesizes(pot1.ctaildom) = pot1.ctailsizes; +nodesizes(pot1.cheaddom) = pot1.cheadsizes; +nodesizes(pot1.ddom) = pot1.dsizes; + +dom_u = mysetdiff(pot2.ctaildom, pot1.cheaddom); +if ~isempty(dom_u) & ~mysubset(dom_u, pot1.ctaildom) + pot1 = extension_pot(pot1, [], [], dom_u, nodesizes(dom_u)); +end + +dom_u = myunion(pot1.cheaddom, pot1.ctaildom); +if ~isempty(dom_u) & ~mysubset(dom_u, pot2.ctaildom) + pot2 = extension_pot(pot2, [], [], dom_u, nodesizes(dom_u)); +end + + +cheaddom = myunion(pot1.cheaddom, pot2.cheaddom); +ctaildom = mysetdiff(myunion(pot1.ctaildom, pot2.ctaildom), cheaddom); +cdom = myunion(cheaddom, ctaildom); +ddom = mysetdiff(domain, cdom); +dsizes = nodesizes(ddom); +dsize = prod(nodesizes(ddom)); +cheadsizes = nodesizes(cheaddom); +cheadsize = sum(nodesizes(cheaddom)); +ctailsizes = nodesizes(ctaildom); +ctailsize = sum(nodesizes(ctaildom)); + +r1 = pot1.cheadsize; +s1 = pot1.ctailsize; +scpot = cell(1, dsize); +mask1 = []; +mask2 = []; +if ~isempty(pot1.ddom) + mask1 = find_equiv_posns(pot1.ddom, ddom); +end +if ~isempty(pot2.ddom) + mask2 = find_equiv_posns(pot2.ddom, ddom); +end +cmask1 = []; +cmask2 = []; +if ~isempty(pot1.cheaddom) + cmask1 = find_equiv_posns(pot1.cheaddom, cheaddom); +end +if ~isempty(pot2.cheaddom) + cmask2 = find_equiv_posns(pot2.cheaddom, cheaddom); +end + +u1 = block(cmask1, cheadsizes); +u2 = block(cmask2, cheadsizes); + +fmaskh = find_equiv_posns(pot1.cheaddom, pot2.ctaildom); +fmaskt = find_equiv_posns(pot1.ctaildom, pot2.ctaildom); + +fh = block(fmaskh, pot2.ctailsizes); +ft = block(fmaskt, pot2.ctailsizes); + +for i=1:dsize + sub = ind2subv(dsizes, i); + sub1 = sub(mask1); + sub2 = sub(mask2); + ind1 = subv2ind(pot1.dsizes, sub1); + ind2 = subv2ind(pot2.dsizes, sub2); + + if isempty(ind1) + ind1 = 1; + end + if isempty(ind2) + ind2 = 1; + end + potc1 = struct(pot1.scgpotc{ind1}); + potc2 = struct(pot2.scgpotc{ind2}); + p = potc1.p; + q = potc2.p; + ro = p*q; + + A = potc1.A; + B = potc1.B; + C = potc1.C; + + E = potc2.A; + F = potc2.B; + G = potc2.C; + + F1 = F(:, fh); + F2 = F(:, ft); + + if ~isempty(F1) + K1 = F1*A; + K2 = F1*B; + FCF = F1*C*F1'; + FC = F1*C; + CFT = C*F1'; + else + K1 = zeros(size(E)); + K2 = zeros(size(F2)); + FCF = zeros(size(G)); + FC = zeros(size(C, 1), size(G, 2)); + CFT = zeros(size(G, 2), size(C, 1)); + end + + + U = zeros(cheadsize,1); + W = zeros(cheadsize,cheadsize); + V = zeros(cheadsize,ctailsize); + + if cheadsize > 0 + U(u1) = A; + U(u2) = E + K1; + W(u1, u1) = C; + W(u2, u2) = G + FCF; + W(u1, u2) = CFT; + W(u2, u1) = FC; + else + U = zeros(cheadsize,1); + W = zeros(cheadsize,cheadsize); + end + if cheadsize > 0 | ctailsize > 0 + if ~isempty(u1) + V(u1, :) = B; + else + V(u1, :) = zeros(potc1.cheadsize, ctailsize); + end + if ~isempty(u2) + V(u2, :) = F2 + K2; + else + V(u2, :) = zeros(potc2.cheadsize, ctailsize); + end + else + V = zeros(cheadsize,ctailsize); + end + + scpot{i} = scgcpot(cheadsize, ctailsize, ro, U, V, W); +end + +pot = scgpot(ddom, cheaddom, ctaildom, nodesizes, scpot); diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/extension_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/extension_pot.m new file mode 100644 index 00000000..a4c2b0e6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/extension_pot.m @@ -0,0 +1,55 @@ +function pot = extension_pot(oldpot, ddom_u, dsizes, ctaildom_u, csizes) +% EXTENSION_POT Extense a stable CG potential. +% pot = extension_pot(oldpot, ddom_u, ctaildom_u, dsizes, csizes) +% ddom_u Added discrete nodes +% ctaildom_u Added continuous tail nodes +% csizes is the size of the tail nodes. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% A CG potential can be extended by adding discrete variables to its % +% domain of continuous variables to its tail % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +ddom = myunion(oldpot.ddom, ddom_u); +ctaildom = myunion(oldpot.ctaildom, ctaildom_u); +cheaddom = oldpot.cheaddom; +udom = myunion(ddom_u, ctaildom_u); +domain = myunion(oldpot.domain, udom); + +ns = zeros(1,max(domain)); +ns(ddom_u) = dsizes; +ns(ctaildom_u) = csizes; +ns(oldpot.ddom) = oldpot.dsizes; +ns(oldpot.cheaddom) = oldpot.cheadsizes; +ns(oldpot.ctaildom) = oldpot.ctailsizes; + +dsizes = ns(ddom); +dsize = prod(ns(ddom)); +cheadsizes = ns(cheaddom); +cheadsize = sum(ns(cheaddom)); +ctailsizes = ns(ctaildom); +ctailsize = sum(ns(ctaildom)); + +BZ = zeros(cheadsize, ctailsize); +potarray = cell(1, dsize); +mask = find_equiv_posns(oldpot.ddom, ddom); + +tmask = find_equiv_posns(oldpot.ctaildom, ctaildom); +tu = block(tmask, ctailsizes); + +for i=1:dsize + sub1 = ind2subv(dsizes, i); + sub2 = sub1(mask); + ind = subv2ind(oldpot.dsizes, sub2); + if isempty(ind) + ind = 1; + end + potc = struct(oldpot.scgpotc{ind}); + p = potc.p; + B = BZ; + if ~isempty(B) + B(:, tu) = potc.B; + end + potarray{i} = scgcpot(cheadsize, ctailsize, p, potc.A, B, potc.C); +end + +pot = scgpot(ddom, cheaddom, ctaildom, ns,potarray); diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/marginalize_pot.m new file mode 100644 index 00000000..317f9f10 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/marginalize_pot.m @@ -0,0 +1,92 @@ +function smallpot = marginalize_pot(bigpot, keep) +% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, keep) + +sumover = mysetdiff(bigpot.domain, keep); +cdom = myunion(bigpot.cheaddom, bigpot.ctaildom); +csumover = myintersect(sumover, bigpot.cheaddom); +dsumover = myintersect(sumover, bigpot.ddom); + +dkeep = myintersect(keep, bigpot.ddom); +ckeep = myintersect(keep, bigpot.cheaddom); +cheaddom = myintersect(keep, bigpot.cheaddom); + +assert(isempty(myintersect(csumover,bigpot.ctaildom))); +ns = zeros(1, max(bigpot.domain)); +ns(bigpot.ddom) = bigpot.dsizes; +ns(bigpot.cheaddom) = bigpot.cheadsizes; +ns(bigpot.ctaildom) = bigpot.ctailsizes; + + +if sum(ns(csumover)) > 0 + for i=1:bigpot.dsize + bigpot.scgpotc{i} = marginalize_pot(bigpot.scgpotc{i}, ckeep, csumover, ns); + end +end + +if (isequal(csumover, cheaddom)) + bigpot.ctaildom = []; +end +% If we are not marginalizing over any discrete nodes, we are done. +if prod(ns(dsumover))==1 + smallpot = scgpot(dkeep, cheaddom, bigpot.ctaildom, ns, bigpot.scgpotc); + return; +end + +if (~isempty(bigpot.ctaildom)) + assert(0); + return; +end + +I = prod(ns(dkeep)); +J = prod(ns(dsumover)); +C = sum(ns(ckeep)); +sum_map = find_equiv_posns(dsumover, bigpot.ddom); +keep_map = find_equiv_posns(dkeep, bigpot.ddom); +iv = zeros(1, length(bigpot.ddom)); % index vector + +p1 = zeros(I,J); +A1 = zeros(C,J,I); +C1 = zeros(C,C,J,I); +for i=1:I + keep_iv = ind2subv(ns(dkeep), i); + iv(keep_map) = keep_iv; + for j=1:J + sum_iv = ind2subv(ns(dsumover), j); + iv(sum_map) = sum_iv; + k = subv2ind(ns(bigpot.ddom), iv); + pot = struct(bigpot.scgpotc{k}); % violate object privacy + p1(i,j) = pot.p; + if C > 0 % so mu1 and Sigma1 are non-empty + A1(:,j,i) = pot.A; + C1(:,:,j,i) = pot.C; + end + end +end + +% Collapse the mixture of Gaussians +coef = mk_stochastic(p1); % coef must be convex combination +%keyboard +p2 = sum(p1,2); +if (all(p2 == 0)) + p2 = p2 + (p2==0)*eps; +end +A = []; +S = []; + +pot = cell(1,I); +ctailsize = sum(ns(bigpot.ctaildom)); +tB = zeros(C, ctailsize); +for i=1:I + if C > 0 + [A, S] = collapse_mog(A1(:,:,i), C1(:,:,:,i), coef(i,:)); + end + p = p2(i); + pot{i} = scgcpot(C, ctailsize, p, A, tB, S); +end + +smallpot = scgpot(dkeep, ckeep, bigpot.ctaildom, ns, pot); + + + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/normalize_pot.m new file mode 100644 index 00000000..da265ef5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/normalize_pot.m @@ -0,0 +1,16 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the SCG potential Pr(X,E) into Pr(X|E) and return log Pr(E). +% [pot, loglik] = normalize_pot(pot) + +% Marginalize down to [], so that the normalizing constant becomes Pr(E) +temp = marginalize_pot(pot, []); +[temp2, loglik] = normalize_pot(temp.scgpotc{1}); + +% Adjust scale factor to reflect the fact that the pot now represents Pr(X | E) instead of Pr(X,E). + +scale = -loglik; +if 1 + for i=1:pot.dsize + pot.scgpotc{i} = rescale_pot( pot.scgpotc{i}, scale); + end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/pot_to_marginal.m new file mode 100644 index 00000000..f7478ef8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/pot_to_marginal.m @@ -0,0 +1,29 @@ +function m = pot_to_marginal(pot) +% POT_TO_MARGINAL Convert a scgpot to a marginal structure. +% m = pot_to_marginal(pot) + +assert(isempty(pot.ctaildom)) +m.domain = pot.domain; +n = pot.cheadsize; +d = pot.dsize; + +if n==0 + m.mu = []; + m.Sigma = []; +else + m.mu = zeros(n, d); + m.Sigma = zeros(n, n, d); +end +%m.T = 0*myones(pot.dsizes); +m.T = 0*myones(pot.dsize); +for i=1:pot.dsize + potc = struct(pot.scgpotc{i}); % violate privacy of object + if n > 0 + m.mu(:,i) = potc.A; + m.Sigma(:,:,i) = potc.C; + end + m.T(i) = potc.p; +end +if isvectorBNT(m.T) + m.T = m.T(:)'; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/recursive_combine_pots.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/recursive_combine_pots.m new file mode 100644 index 00000000..27c04b75 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/recursive_combine_pots.m @@ -0,0 +1,40 @@ +function pot = recursive_combine_pots(pot1, pot2) +% RECURSIVE_COMBINE_POTS recursive combine two potentials +% pot = recursive_combine_pots(pot1, pot2) + +pot1 = reduce_pot(pot1); +pot2 = reduce_pot(pot2); +% Recursion is stopped, if recusive-combination is defined by direct combination, +% i.e. if the domain of one potential is disjoint from the head of the other. +if (isempty(myintersect(pot1.domain,pot2.cheaddom))|... + isempty(myintersect(pot1.cheaddom,pot2.domain))) + pot = direct_combine_pots(pot1,pot2); +else + % Test wether one of the set-differences is not empty + % as defined in Lauritzen99 "Stable Local Computation with Conditional Gaussian Distributions" + % on page 9 + D12 = mysetdiff(pot1.cheaddom, pot2.domain); + D21 = mysetdiff(pot2.cheaddom, pot1.domain); + if (isempty(D12) & isempty(D21)) + assert(0,'Recursive combination is not defined'); + end + + if ~isempty(D12) + % Calculate the complementary potential for the set + % D1\D12 as defined in Lauritzen 99, page 9 + keep = mysetdiff(pot1.domain,D12); + [margpot, comppot] = complement_pot(pot1,keep); + margpot = reduce_pot(margpot); + comppot = reduce_pot(comppot); + pot = direct_combine_pots( recursive_combine_pots(margpot, pot2), comppot); + elseif ~isempty(D21) + keep = mysetdiff(pot2.domain,D21); + [margpot, comppot] = complement_pot(pot2,D21); + margpot = reduce_pot(margpot); + comppot = reduce_pot(comppot); + pot = direct_combine_pots( recursive_combine_pots(pot1, margpot), comppot); + end +end + + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/reduce_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/reduce_pot.m new file mode 100644 index 00000000..7ee4ab49 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/reduce_pot.m @@ -0,0 +1,53 @@ +function [reduced_pot,successful] = reduce_pot(pot,tailnodes) +% Executes the reduce operation defined in +% Stable Local Computation with Conditional Gaussian Distributions +% Steffen L. Lauritzen +% Frank Jensen +% September 1999 +% The potential pot is reduced if B contains any zero columns +% The test are restricted to the positions in tailnodes. +% Any columns successfully deleted are entered in the array successful +if nargin < 2 + tailnodes = pot.ctaildom; +end + +successful = []; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Keep track of remaining tailnodes % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +rem_tailnodes = pot.ctaildom; +for i = tailnodes + pos = find(i==rem_tailnodes); + successful_red = [pos]; + red_scgcpot = cell(1,pot.dsize); + j = 1; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Test whether all components of pot.scgpotc can be reduced % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + while ((j <= pot.dsize) & ~isempty(successful_red)) + [cpot,successful_red] = reduce_pot(pot.scgpotc{j},pos); + red_scgcpot{j} = cpot; + j = j + 1; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % If i is a reducible tailnode, then reduce the potential % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if ~isempty(successful_red) + successful = [successful i]; + pot.scgpotc = red_scgcpot; + rem_tailnodes = mysetdiff(rem_tailnodes,i); + end; +end + +pot.ctaildom = rem_tailnodes; +positions = find_equiv_posns(rem_tailnodes,pot.ctaildom); +pot.ctailsizes = pot.ctailsizes(positions); +pot.ctailsize = sum(pot.ctailsizes); +pot.domain = mysetdiff(pot.domain,successful); +reduced_pot = pot; + + + + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/scgpot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/scgpot.m new file mode 100644 index 00000000..7f923f22 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/scgpot.m @@ -0,0 +1,30 @@ +function pot = scgpot(ddom, cheaddom, ctaildom, node_sizes, scgpotc) +% SCGPOT Make a stable CG potential. +% pot = scgpot(ddom, cheaddom, ctaildom, node_sizes, scgpotc) +% +% ddom is discrete nodes contains in the potential +% cheaddom is head nodes constains in the potential +% ctaildom is tail nodes contains in the potential +% node_sizes(i) is the size of the i'th node. +% scgpotc is list of scgcpot objects. + +pot.ddom = ddom; +pot.cheaddom = cheaddom; +pot.ctaildom = ctaildom; +pot.domain = myunion(ddom, myunion(cheaddom, ctaildom)); +pot.dsizes = node_sizes(pot.ddom); +pot.dsize = prod(node_sizes(pot.ddom)); +pot.cheadsizes = node_sizes(pot.cheaddom); +pot.cheadsize = sum(node_sizes(pot.cheaddom)); +pot.ctailsizes = node_sizes(pot.ctaildom); +pot.ctailsize = sum(node_sizes(pot.ctaildom)); + +if nargin < 5 + scgpotc = cell(1, pot.dsize); + for i=1:pot.dsize + scgpotc{i} = scgcpot(pot.cheadsize, pot.ctailsize); + end +end +pot.scgpotc = scgpotc; + +pot = class(pot, 'scgpot'); diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Entries new file mode 100644 index 00000000..25bc0a31 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Entries @@ -0,0 +1,10 @@ +/approxeq_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/display.m/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/multiply_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002// +/upot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/upot_to_opt_policy.m/1.1.1.1/Wed May 29 15:59:58 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Repository new file mode 100644 index 00000000..b7cf8acd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/@upot diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/approxeq_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/approxeq_pot.m new file mode 100644 index 00000000..de47ee12 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/approxeq_pot.m @@ -0,0 +1,5 @@ +function p = approxeq_pot(A, B, tol) + +if nargin < 3, tol = 1e-3; end + +p = approxeq(A.p, B.p, tol) & approxeq(A.u, B.u, tol); diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/display.m b/sourcecodes/bnt-master/BNT/potentials/@upot/display.m new file mode 100644 index 00000000..bb8e76e9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/display.m @@ -0,0 +1,4 @@ +function display(pot) + +disp('utility potential object'); +disp(struct(pot)); diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/divide_by_pot.m new file mode 100644 index 00000000..8ac75785 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/divide_by_pot.m @@ -0,0 +1,13 @@ +function Tbig = divide_by_pot(Tbig, Tsmall) +% DIVIDE_BY_POT Tbig /= Tsmall +% Tbig = divide_by_pot(Tbig, Tsmall) +% +% Tsmall's domain must be a subset of Tbig's domain. + +smallp = extend_domain_table(Tsmall.p, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes); +smallp = smallp + (smallp==0); +Tbig.p = Tbig.p ./ smallp; + +smallu = extend_domain_table(Tsmall.u, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes); +Tbig.u = Tbig.u - smallu; + diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/marginalize_pot.m new file mode 100644 index 00000000..be4bf249 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/marginalize_pot.m @@ -0,0 +1,18 @@ +function smallpot = marginalize_pot(bigpot, onto, maximize) +% MARGINALIZE_POT Marginalize a upot onto a smaller domain. +% smallpot = marginalize_pot(bigpot, onto, maximize) +% +% The maximize argument is ignored + +numer = marg_table(bigpot.p .* bigpot.u, bigpot.domain, bigpot.sizes, onto); +denom = marg_table(bigpot.p, bigpot.domain, bigpot.sizes, onto); + +p = denom; +% replace 0s by 1s before dividing. This is valid since demon(i) = 0 => numer(i) = 0 +denom = denom + (denom == 0); +u = numer ./ denom; + +ns = zeros(1, max(bigpot.domain)); +ns(bigpot.domain) = bigpot.sizes; + +smallpot = upot(onto, ns(onto), p, u); diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/multiply_by_pot.m new file mode 100644 index 00000000..28b3586d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/multiply_by_pot.m @@ -0,0 +1,12 @@ +function Tbig = multiply_by_pot(Tbig, Tsmall) +% MULTIPLY_BY_POT Tbig *= Tsmall +% Tbig = multiply_by_pot(Tbig, Tsmall) +% +% Tsmall's domain must be a subset of Tbig's domain. + +smallp = extend_domain_table(Tsmall.p, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes); +Tbig.p = Tbig.p .* smallp; + +smallu = extend_domain_table(Tsmall.u, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes); +Tbig.u = Tbig.u + smallu; + diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/normalize_pot.m new file mode 100644 index 00000000..dc3dfa8a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/normalize_pot.m @@ -0,0 +1,13 @@ +function [pot, loglik] = normalize_pot(pot) +% NORMALIZE_POT Convert the probability part of a utility potential +% [pot, loglik] = normalize_pot(pot) + +[pot.p, lik] = normalise(pot.p); +%pot.u = pot.u - sum(pot.u(:)); +%pot.u = pot.u ./ sum(pot.u(:)); % same as normalise(pot.u) +%pot.u = normalise(pot.u); +%pot.u = pot.u / 726.8121; +pot.u = pot.u / 10; +loglik = log(lik + (lik==0)*eps); + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@upot/pot_to_marginal.m new file mode 100644 index 00000000..4d3a1c37 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/pot_to_marginal.m @@ -0,0 +1,9 @@ +function m = pot_to_marginal(pot) +% POT_TO_MARGINAL Convert a upot to a structure. +% m = pot_to_marginal(pot) + +m.domain = pot.domain; +m.T = pot.p; +m.U = pot.u; + + diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/upot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/upot.m new file mode 100644 index 00000000..502858b8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/upot.m @@ -0,0 +1,15 @@ +function pot = upot(domain, sizes, p, u) +% UPOT Make a discrete utility potential. +% pot = upot(domain, sizes, p, u) +% +% sizes(i) is the size of the i'th domain element. +% p defaults to all 1s, u defaults to all 0s. + +if nargin < 3, p = myones(sizes); end +if nargin < 4, u = 0*myones(sizes); end + +pot.domain = domain; +pot.p = myreshape(p, sizes); +pot.u = myreshape(u, sizes); +pot.sizes = sizes(:)'; +pot = class(pot, 'upot'); diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/upot_to_opt_policy.m b/sourcecodes/bnt-master/BNT/potentials/@upot/upot_to_opt_policy.m new file mode 100644 index 00000000..20f8d2ec --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/@upot/upot_to_opt_policy.m @@ -0,0 +1,25 @@ +function [policy, EU] = upot_to_opt_policy(pot) +% UPOT_TO_OPT_POLICY Compute an optimal deterministic policy given a utility potential +% [policy, EU] = upot_to_opt_policy(pot) +% +% policy(a,b, ..., z) = P(do z | a, b, ..), which will be a delta function +% EU is the contraction of this potential, i.e., P .* U + +sz = pot.sizes; % mysize(pot.p); +if isempty(sz) + EU = pot.u; + policy = []; + return; +end + +parent_size = prod(sz(1:end-1)); +self_size = sz(end); +C = pot.p .* pot.u; % contraction +C = reshape(C, parent_size, self_size); +policy = zeros(parent_size, self_size); +for i=1:parent_size + act = argmax(C(i,:)); + policy(i, act) = 1; +end +policy = myreshape(policy, sz); +EU = sum(C(:)); diff --git a/sourcecodes/bnt-master/BNT/potentials/CPD_to_pot.m b/sourcecodes/bnt-master/BNT/potentials/CPD_to_pot.m new file mode 100644 index 00000000..8d0aa3ce --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/CPD_to_pot.m @@ -0,0 +1,21 @@ +function pot = CPD_to_pot(pot_type, CPD, domain, ns, cnodes, evidence) +% CPD_TO_POT Convert a CPD to a potential of the specified form, incorporating any evidence +% pot = CPD_to_pot(pot_type, CPD, domain, node_sizes, cnodes, evidence) +% +% pot_type is one of 'd', 'g', or 'cg'. +% domain is the domain of CPD. +% node_sizes(i) is the size of node i. +% cnodes = the cts nodes +% evidence{i} is the evidence on the i'th node. + +switch pot_type + case 'd', + pot = CPD_to_dpot(CPD, domain, ns, cnodes, evidence); + case 'g', + pot = CPD_to_cpot(CPD, domain, ns, cnodes, evidence); + case 'cg', + pot = CPD_to_cgpot(CPD, domain, ns, cnodes, evidence); + otherwise, + error(['can''t handle pot_type ' pot_type]); +end + diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries new file mode 100644 index 00000000..bfef35ea --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries @@ -0,0 +1,7 @@ +/CPD_to_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +/README/1.1.1.1/Wed May 29 15:59:58 2002// +/check_for_cd_arcs.m/1.1.1.1/Wed May 29 15:59:58 2002// +/determine_pot_type.m/1.1.1.1/Wed May 29 15:59:58 2002// +/genops.c/1.1.1.1/Sat Jul 28 15:43:40 2001// +/mk_initial_pot.m/1.1.1.1/Wed May 29 15:59:58 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries.Log new file mode 100644 index 00000000..35415310 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries.Log @@ -0,0 +1,9 @@ +A D/@cgpot//// +A D/@cpot//// +A D/@dpot//// +A D/@mpot//// +A D/@scgcpot//// +A D/@scgpot//// +A D/@upot//// +A D/Old//// +A D/Tables//// diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/CVS/Repository new file mode 100644 index 00000000..3517a1ca --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Entries new file mode 100644 index 00000000..e7d041f1 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Entries @@ -0,0 +1,2 @@ +/comp_eff_node_sizes.m/1.1.1.1/Wed May 29 15:59:58 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Repository new file mode 100644 index 00000000..48670656 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/Old diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/comp_eff_node_sizes.m b/sourcecodes/bnt-master/BNT/potentials/Old/comp_eff_node_sizes.m new file mode 100644 index 00000000..e245f602 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Old/comp_eff_node_sizes.m @@ -0,0 +1,11 @@ +function ens = comp_eff_node_sizes(ns, cnodes, ev, domain) + +dnodes = mysetdiff(1:length(ns), cnodes); +odom = domain(~isemptycell(evidence(domain))); +cdom = myintersect(cnodes, domain); +ddom = myintersect(dnodes, domain); +cobs = myintersect(cdom, odom); +dobs = myintersect(ddom, odom); +ens = ns; +ens(cobs) = 0; +ens(dobs) = 1; diff --git a/sourcecodes/bnt-master/BNT/potentials/README b/sourcecodes/bnt-master/BNT/potentials/README new file mode 100644 index 00000000..116266ee --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/README @@ -0,0 +1,100 @@ +The following kinds of potentials are supported +- dpot: discrete +- upot: utility +- mpot: Gaussian in moment form +- cpot: Gaussian in canonical form +- cgpot: conditional (mixture) Gaussian, a list of mpots/cpot +- scgpot: stable conditional Gaussian, a list of scgcpots +- scgcpot: just used by scgpot + +Many of these are described in the following book + +@book{Cowell99, + author = "R. G. Cowell and A. P. Dawid and S. L. Lauritzen and D. J. Spiegelhalter", + title = "Probabilistic Networks and Expert Systems", + year = 1999, + publisher = "Springer" +} + +CPD_to_pot converts P(Z|A,B,...) to phi(A,B,...,Z). + +A table is like a dpot, except it is a structure, not an object. +Code that uses tables is faster but less flexible. + + ----------- + +A potential is a joint probability distribution on a set of nodes, +which we call the potential's domain (which is always sorted). +A potential supports the operations of multiplication and +marginalization. + +If the nodes are discrete, the potential can be represented as a table +(multi-dimensional array). If the nodes are Gaussian, the potential +can be represented as a quadratic form. If there are both discrete and +Gaussian nodes, we use a table of quadratic forms. For details on the +Gaussian case, see below. + +For discrete potentials, the 'sizes' field specifies the number of +values each node in the domain can take on. For continuous potentials, +the 'sizes' field specifies the block-size of each node. + +If some of the nodes are observed, extra complications arise. We +handle the discrete and continuous cases differently. Suppose the +domain is [X Y], with sizes [6 2], where X is observed to have value x. +In the discrete case, the potential will have many zeros in it +(T(X,:) will be 0 for all X ~= x), which can be inefficient. Instead, +we set sizes to [1 2], to indicate that X has only one possible value +(namely x). For continuous nodes, we set sizes = [0 2], to indicate that X no +longer appears in the mean vector or covariance matrix (we must avoid +0s in Sigma, lest it be uninvertible). When a potential is created, we +assume the sizes of the nodes have been adjusted to include the +evidence. This is so that the evidence can be incorporated at the +outset, and thereafter the inference algorithms can ignore it. + + ------------ + +A Gaussian potential can be represented in terms of its +moment characteristics (mu, Sigma, logp), or in terms of its canonical +characteristics (g, h, K). Although the moment characteristics are +more familiar, it turns out that canonical characteristics are +more convenient for the junction tree algorithm, for the same kinds of +reasons why backwards inference in an LDS uses the information form of +the Kalman filter (see Murphy (1998a) for a discussion). + +When working with *conditional* Gaussian potentials, the method proposed +by Lauritzen (1992), and implemented here, requires converting from +canonical to moment form before marginalizing the discrete variables, +and converting back from moment to canonical form before +multiplying/dividing. A new algorithm, due to Lauritzen and Jensen +(1999), works exclusively in moment form, and +hence is more numerically stable. It can also handle 0s in the +covariance matrix, i.e., deterministic relationships between cts +variables. However, it has not yet been implemented, +since it requires major changes to the jtree algorithm. + +In Murphy (1998b) we extend Lauritzen (1992) to handle +vector-valued nodes. This means the vectors and matrices become block +vectors and matrices. This manifests itself in the code as in the +following example. +Suppose we have a potential on nodes dom=[3,4,7] with block sizes=[2,1,3]. +Then nodes 3 and 7 correspond to blocks 1,3 which correspond to indices 1,2,4,5,6. +>> find_equiv_posns([3 7], dom)=[1,3] +>> block([1,3],blocks)=[1,2,4,5,6]. + +For more details, see + +- "Filtering and Smoothing in Linear Dynamical Systems using the Junction Tree Algorithm", + K. Murphy, 1998a. UCB Tech Report. + +- "Inference and learning in hybrid Bayesian networks", + K. Murphy. UCB Technical Report CSD-98-990, 1998b. + +- "Propagation of probabilities, means and variances in mixed + graphical association models", S. L. Lauritzen, 1992, JASA 87(420):1098--1108. + +- "Causal probabilistic networks with both discrete and continuous variables", + K. G. Olesen, 1993. PAMI 3(15). This discusses implementation details. + +- "Stable local computation with Conditional Gaussian distributions", + S. Lauritzen and F. Jensen, 1999. Univ. Aalborg Tech Report R-99-2014. + www.math.auc.dk/research/Reports.html. diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries new file mode 100644 index 00000000..f550f801 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries @@ -0,0 +1,19 @@ +/divide_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_table.m/1.1.1.1/Thu Aug 5 15:25:54 2004// +/extend_domain_table.m/1.1.1.1/Wed Aug 4 15:53:26 2004// +/marg_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/marg_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/marg_table.m/1.1.1.1/Wed Aug 4 15:51:48 2004// +/marg_tableC.c/1.1.1.1/Wed Oct 2 15:39:02 2002// +/marg_tableM.m/1.1.1.1/Tue Oct 1 17:39:08 2002// +/mult_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/mult_by_table.c/1.1.1.1/Tue Oct 1 21:23:22 2002// +/mult_by_table.m/1.1.1.1/Wed Aug 4 15:52:58 2004// +/mult_by_table2.m/1.1.1.1/Wed Oct 2 15:30:32 2002// +/mult_by_tableC.c/1.1.1.1/Tue Oct 1 21:33:50 2002// +/mult_by_tableM.m/1.1.1.1/Wed Oct 2 15:28:48 2002// +/mult_by_table_global.m/1.1.1.1/Sun Sep 29 10:21:30 2002// +/rep_mult.c/1.1.1.1/Wed May 29 15:59:58 2002// +/repmat_and_mult.c/1.1.1.1/Tue Oct 1 21:20:00 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository new file mode 100644 index 00000000..e52428a3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/Tables diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c new file mode 100644 index 00000000..d7645ab1 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c @@ -0,0 +1,126 @@ +/* divide_by_sparse_table.c ../potential/tables*/ + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex; + int *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + + plhs[0] = mxDuplicateArray(prhs[0]); + pbDomain = mxGetPr(prhs[1]); + bdim = mxGetNumberOfElements(prhs[1]); + psDomain = mxGetPr(prhs[4]); + sdim = mxGetNumberOfElements(prhs[4]); + + pbSize = mxGetPr(prhs[2]); + psSize = mxGetPr(prhs[5]); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + bpr = mxGetPr(plhs[0]); + bir = mxGetIr(plhs[0]); + bjc = mxGetJc(plhs[0]); + NZB = bjc[1]; + + spr = mxGetPr(prhs[3]); + sir = mxGetIr(prhs[3]); + sjc = mxGetJc(prhs[3]); + NZS = sjc[1]; + + if(sdim == 0){ + value = *spr; + if(value == 0)value = 1; + for(i=0; i<NZB; i++){ + bpr[i] /= value; + } + return; + } + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + bpr[i] /= spr[position]; + } + } + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c new file mode 100644 index 00000000..6942eddd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c @@ -0,0 +1,120 @@ +/* divide_by_table.c ../potential/tables */ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2; + double *pbDomain, *psDomain, *sp, *zp, *bs, value; + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[4]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[4]); + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[3]); + sp = mxGetPr(prhs[3]); + + bs = mxGetPr(prhs[2]); + + if(NS == 1){ + value = *sp; + if(value == 0) value = 1; + for(i=0; i<NB; i++){ + zp[i] /= value; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + value = sp[i]; + if(value == 0) value = 1; + zp[i] /= value; + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)bs[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + value = *sp; + if(value == 0) value = 1; + *zp++ /= value; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m new file mode 100644 index 00000000..6f10718b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m @@ -0,0 +1,12 @@ +function bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% DIVIDE_BY_TABLE +% bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +% Replace 0s by 1s before dividing. This is valid, Ts(i)=0 iff Tbig(i)=0. +Ts = Ts + (Ts==0); +%Tbig.T(:) = Tbig.T(:) ./ Ts(:); +bigT(:) = bigT(:) ./ Ts(:); + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m new file mode 100644 index 00000000..4ed0b2cd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m @@ -0,0 +1,24 @@ +function B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz) +% EXTEND_DOMAIN_TABLE Expand an array so it has the desired size. +% B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz) +% +% A is the array with domain smalldom and sizes smallsz. +% bigdom is the desired domain, with sizes bigsz. +% +% Example: +% smalldom = [1 3], smallsz = [2 4], bigdom = [1 2 3 4], bigsz = [2 1 4 5], +% so B(i,j,k,l) = A(i,k) for i in 1:2, j in 1:1, k in 1:4, l in 1:5 + +if isequal(size(A), [1 1]) % a scalar + B = A; % * myones(bigsz); + return; +end + +map = find_equiv_posns(smalldom, bigdom); +sz = ones(1, length(bigdom)); +sz(map) = smallsz; +B = myreshape(A, sz); % add dimensions for the stuff not in A +sz = bigsz; +sz(map) = 1; % don't replicate along A's dimensions +B = myrepmat(B, sz(:)'); + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c new file mode 100644 index 00000000..f6f7aa26 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c @@ -0,0 +1,181 @@ +/* marg_sparse_table.c ../potential/tables*/ + +/******************************************/ +/* 5 input & 1 output */ +/* Big sparse table */ +/* Big domain */ +/* Big sizes */ +/* onto */ +/* maximize, if missed, maximize=0 */ +/* */ +/* small sparse table */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +mxArray* convert_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = Table[temp]; + count++; + } + return spTable; +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, bdim, sdim, NS, NZB, position, bindex, sindex, maximize, nzCounts=0; + int *mask, *sequence, *result, *bir, *bjc, *ssize, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *sTable, *pbDomain, *psDomain, *pbSize, *bpr, *spr; + const char *field_names[] = {"domain", "T", "sizes"}; + + if(nrhs < 5) maximize = 0; + else maximize = (int)mxGetScalar(prhs[4]); + + bdim = mxGetNumberOfElements(prhs[1]); + sdim = mxGetNumberOfElements(prhs[3]); + pbSize = mxGetPr(prhs[2]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[3]); + bpr = mxGetPr(prhs[0]); + bir = mxGetIr(prhs[0]); + bjc = mxGetJc(prhs[0]); + NZB = bjc[1]; + + if(sdim == 0){ + plhs[0] = mxCreateSparse(1, 1, 1, mxREAL); + spr = mxGetPr(plhs[0]); + bir = mxGetIr(plhs[0]); + bjc = mxGetJc(plhs[0]); + *spr = 0; + *bir = 0; + bjc[0] = 0; + bjc[1] = 1; + if(maximize){ + for(i=0; i<NZB; i++){ + *spr = (*spr < bpr[i])? bpr[i] : *spr; + } + } + else{ + for(i=0; i<NZB; i++){ + *spr += bpr[i]; + } + } + return; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + sTable = malloc(NZB * sizeof(double)); + sequence = malloc(NZB * 2 * sizeof(double)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + ssize = malloc(sdim * sizeof(int)); + + NS = 1; + for(i=0; i<count; i++){ + ssize[i] = (int)pbSize[mask[i]]; + NS *= ssize[i]; + } + + for(i=0; i<NZB; i++)sTable[i] = 0; + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * ssize[i]; + } + + count = 0; + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare); + if(result){ + position = (result - sequence) / 2; + if(maximize) + sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position]; + else sTable[position] += bpr[i]; + } + else { + if(maximize) + sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts]; + else sTable[nzCounts] += bpr[i]; + sequence[count] = sindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + plhs[0] = convert_table_to_sparse(sTable, sequence, nzCounts, NS); + + free(sTable); + free(sequence); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); + free(ssize); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c new file mode 100644 index 00000000..16329dee --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c @@ -0,0 +1,175 @@ +/* marg_table.c ../potential/tables */ + + +/******************************************/ +/* 5 input & 1 output */ +/* Big table */ +/* Big domain */ +/* Big sizes */ +/* onto */ +/* maximize, if missed, maximize=0 */ +/* */ +/* small table */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize; + double *pb, *ps, *bp, *sp, *pbd; + + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[3]); + pb = mxGetPr(prhs[1]); + ps = mxGetPr(prhs[3]); + + NB = mxGetNumberOfElements(prhs[0]); + bp = mxGetPr(prhs[0]); + + pbd = mxGetPr(prhs[2]); + + if(nrhs < 5) maximize = 0; + else maximize = (int)mxGetScalar(prhs[4]); + + if(siz_s == 0){ + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); + sp = mxGetPr(plhs[0]); + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + ssize = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(ps[i] == pb[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbd[i]; + sy[i] = 1; + } + for(i=0; i<siz_s; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + ssize[i] = sx[temp]; + } + + NS = 1; + for(i=0; i<ndim; i++){ + NS *= sy[i]; + } + + plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL); + sp = mxGetPr(plhs[0]); + + if(NS == 1){ + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++) *sp++ = *bp++; + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + if(maximize){ + for(j=0; j<NB; j++){ + *sp = (*sp < *bp)? *bp : *sp; + bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + else{ + for(j=0; j<NB; j++){ + *sp += *bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); + free(ssize); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m new file mode 100644 index 00000000..0f1e5be5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m @@ -0,0 +1,27 @@ +function smallT = marg_table(bigT, bigdom, bigsz, onto, maximize) +% MARG_TABLE Marginalize a table +% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize) + +if nargin < 5, maximize = 0; end + + +smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array +sum_over = mysetdiff(bigdom, onto); +ndx = find_equiv_posns(sum_over, bigdom); +if maximize + for i=1:length(ndx) + smallT = max(smallT, [], ndx(i)); + end +else + for i=1:length(ndx) + smallT = sum(smallT, ndx(i)); + end +end + + +ns = zeros(1, max(bigdom)); +%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1 +ns(bigdom) = bigsz; + +smallT = squeeze(smallT); % remove all dimensions of size 1 +smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1 diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c new file mode 100644 index 00000000..16329dee --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c @@ -0,0 +1,175 @@ +/* marg_table.c ../potential/tables */ + + +/******************************************/ +/* 5 input & 1 output */ +/* Big table */ +/* Big domain */ +/* Big sizes */ +/* onto */ +/* maximize, if missed, maximize=0 */ +/* */ +/* small table */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize; + double *pb, *ps, *bp, *sp, *pbd; + + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[3]); + pb = mxGetPr(prhs[1]); + ps = mxGetPr(prhs[3]); + + NB = mxGetNumberOfElements(prhs[0]); + bp = mxGetPr(prhs[0]); + + pbd = mxGetPr(prhs[2]); + + if(nrhs < 5) maximize = 0; + else maximize = (int)mxGetScalar(prhs[4]); + + if(siz_s == 0){ + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); + sp = mxGetPr(plhs[0]); + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + ssize = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(ps[i] == pb[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbd[i]; + sy[i] = 1; + } + for(i=0; i<siz_s; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + ssize[i] = sx[temp]; + } + + NS = 1; + for(i=0; i<ndim; i++){ + NS *= sy[i]; + } + + plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL); + sp = mxGetPr(plhs[0]); + + if(NS == 1){ + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++) *sp++ = *bp++; + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + if(maximize){ + for(j=0; j<NB; j++){ + *sp = (*sp < *bp)? *bp : *sp; + bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + else{ + for(j=0; j<NB; j++){ + *sp += *bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); + free(ssize); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m new file mode 100644 index 00000000..22487dff --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m @@ -0,0 +1,30 @@ +function smallT = marg_tableM(bigT, bigdom, bigsz, onto, maximize) +% MARG_TABLE Marginalize a table +% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize) + + % marg_tableM is the same as marg_table.m, but we are sure + % it is not marg_table.c + +if nargin < 5, maximize = 0; end + + +smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array +sum_over = mysetdiff(bigdom, onto); +ndx = find_equiv_posns(sum_over, bigdom); +if maximize + for i=1:length(ndx) + smallT = max(smallT, [], ndx(i)); + end +else + for i=1:length(ndx) + smallT = sum(smallT, ndx(i)); + end +end + + +ns = zeros(1, max(bigdom)); +%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1 +ns(bigdom) = bigsz; + +smallT = squeeze(smallT); % remove all dimensions of size 1 +smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1 diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c new file mode 100644 index 00000000..c1331508 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c @@ -0,0 +1,155 @@ +/* mult_by_sparse_table.c ../potential/tables*/ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ + double *ptr; + void *newptr; + int *ir, *jc; + int nbytes; + + if(new_nzmax == old_nzmax) return; + nbytes = new_nzmax * sizeof(*ptr); + ptr = mxGetPr(spArray); + newptr = mxRealloc(ptr, nbytes); + mxSetPr(spArray, newptr); + nbytes = new_nzmax * sizeof(*ir); + ir = mxGetIr(spArray); + newptr = mxRealloc(ir, nbytes); + mxSetIr(spArray, newptr); + jc = mxGetJc(spArray); + jc[0] = 0; + jc[1] = new_nzmax; + mxSetNzmax(spArray, new_nzmax); +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr; + + pbDomain = mxGetPr(prhs[1]); + bdim = mxGetNumberOfElements(prhs[1]); + psDomain = mxGetPr(prhs[4]); + sdim = mxGetNumberOfElements(prhs[4]); + + pbSize = mxGetPr(prhs[2]); + psSize = mxGetPr(prhs[5]); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + bpr = mxGetPr(prhs[0]); + bir = mxGetIr(prhs[0]); + bjc = mxGetJc(prhs[0]); + NZB = bjc[1]; + + spr = mxGetPr(prhs[3]); + sir = mxGetIr(prhs[3]); + sjc = mxGetJc(prhs[3]); + NZS = sjc[1]; + + plhs[0] = mxDuplicateArray(prhs[0]); + rpr = mxGetPr(plhs[0]); + rir = mxGetIr(plhs[0]); + rjc = mxGetJc(plhs[0]); + rjc[0] = 0; + rjc[1] = NZB; + + if(sdim == 0){ + for(i=0; i<NZB; i++){ + rpr[i] *= *spr; + } + return; + } + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + rpr[nzCounts] = bpr[i] * spr[position]; + rir[nzCounts] = bindex; + nzCounts++; + } + } + + reset_nzmax(plhs[0], NZB, nzCounts); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c new file mode 100644 index 00000000..3f67983d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c @@ -0,0 +1,114 @@ +/* mult_by_table.c ../potential/tables */ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2; + double *pbDomain, *psDomain, *sp, *zp, *bs; + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[4]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[4]); + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[3]); + sp = mxGetPr(prhs[3]); + + bs = mxGetPr(prhs[2]); + + if(NS == 1){ + for(i=0; i<NB; i++){ + zp[i] *= *sp; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + zp[i] *= sp[i]; + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)bs[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + *zp++ *= *sp; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m new file mode 100644 index 00000000..a92e340f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m @@ -0,0 +1,7 @@ +function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% MULT_BY_TABLE +% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m new file mode 100644 index 00000000..92fb4ac7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m @@ -0,0 +1,30 @@ +function bigT = mult_by_table2(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% MULT_BY_TABLE +% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + +%Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +%bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape + +% extend_domain_table has a lot of overhead for small tables, +% since it calls myreshape and myrepmat, which check for 1 dimensional case. +% Here, we check up front. + +if length(bigdom)==1 % vector + bigT = bigT .* smallT; % smallT can be scalar or vector +else + if (length(bigsz) == length(smallsz)) & all(bigsz == smallsz) + bigT = bigT .* smallT; + else + map = find_equiv_posns(smalldom, bigdom); + sz = ones(1, length(bigdom)); + sz(map) = smallsz; + smallT = reshape(smallT, sz); % add dimensions of size 1 for missing domain + % we can use reshape instead of myreshape, because we know length(sz)>1 + sz = bigsz; + sz(map) = 1; % don't replicate along small domain, which is shared + % we can use repmat instead of myrepmat, because we know length(sz)>1 + smallT = repmat(smallT, sz(:)'); + bigT(:) = bigT(:) .* smallT(:); + end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c new file mode 100644 index 00000000..3f67983d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c @@ -0,0 +1,114 @@ +/* mult_by_table.c ../potential/tables */ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2; + double *pbDomain, *psDomain, *sp, *zp, *bs; + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[4]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[4]); + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[3]); + sp = mxGetPr(prhs[3]); + + bs = mxGetPr(prhs[2]); + + if(NS == 1){ + for(i=0; i<NB; i++){ + zp[i] *= *sp; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + zp[i] *= sp[i]; + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)bs[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + *zp++ *= *sp; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m new file mode 100644 index 00000000..a92e340f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m @@ -0,0 +1,7 @@ +function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% MULT_BY_TABLE +% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m new file mode 100644 index 00000000..cc6d7124 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m @@ -0,0 +1,7 @@ +function mult_by_table_global(bigT, bigdom, bigsz, smallT, smalldom, smallsz) + +% all arguments are read only +global NEWBIGT_GLOBAL + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +NEWBIGT_GLOBAL = bigT(:) .* Ts(:); diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c new file mode 100644 index 00000000..c7b687a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c @@ -0,0 +1,92 @@ +/* rep_mult.c repmat first two operands to the size provided by */ +/* the third operand, then perform point multiply */ +/* 3 input, 1 output */ +/* C = rep_mult(A, B, sizes) */ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + double *xp, *yp, *zp, *pSizes; + int xnd, ynd, numElements = 1; + const int *xdim, *ydim; + int i, j, ndim; + int *s, *sx, *sy, *cpsx, *cpsy; + int *subs, *s1, *cpsx2, *cpsy2; + + if (nrhs != 3) + mexErrMsgTxt("Incorrect number of inputs."); + + if (nlhs > 1) + mexErrMsgTxt("Too many output arguments."); + + xnd = mxGetNumberOfDimensions(prhs[0]); + ynd = mxGetNumberOfDimensions(prhs[1]); + xdim = mxGetDimensions(prhs[0]); + ydim = mxGetDimensions(prhs[1]); + ndim = mxGetNumberOfElements(prhs[2]); + + pSizes = mxGetPr(prhs[2]); + + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + s = (int *)malloc(sizeof(int)*ndim); + s1 = (int *)malloc(sizeof(int)*ndim); + *(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1; + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsx2 = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + subs[i] = 0; + sx[i] = (i < xnd) ? xdim[i] : 1; + sy[i] = (i < ynd) ? ydim[i] : 1; + s[i] = (int)pSizes[i]; + s1[i] = s[i] - 1; + numElements *= s[i]; + } + + for(i=0; i<ndim-1; i++){ + cpsx[i+1] = cpsx[i]*sx[i]--; + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsx2[i] = cpsx[i]*sx[i]; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]); + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + zp = mxGetPr(plhs[0]); + xp = mxGetPr(prhs[0]); + yp = mxGetPr(prhs[1]); + + for(j=0; j<numElements; j++){ + *zp++ = *xp * *yp; + for(i=0; i<ndim; i++){ + if(subs[i] == s1[i]){ + subs[i] = 0; + if(sx[i]) + xp -= cpsx2[i]; + if(sy[i]) + yp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sx[i]) + xp += cpsx[i]; + if(sy[i]) + yp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(s1); + free(cpsx); + free(cpsy); + free(subs); + free(cpsx2); + free(cpsy2); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c new file mode 100644 index 00000000..a7bdba35 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c @@ -0,0 +1,97 @@ +/**************************************************** +A = mult_by_array(big, small) +implicitely copies small |big|/|small| times +and then does element-wise multiplication. + +i.e., +C = repmat(small(:), 1, length(big(:))/length(small(:))) +A = reshape(big(:) .* C(:), size(big)) + +However, this C version avoids the expense of the repmat. + +Written by wei.hu@intel.com, 28 Jan 2002. +/****************************************************/ + + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + double *sp, *zp; + int i, j, NB, NS, xnd, ynd, ndim; + const int *xdim, *ydim; + int *s, *sx, *sy, *cpsy, *subs, *cpsy2; + + if (nrhs != 2) + mexErrMsgTxt("Incorrect number of inputs."); + + if (nlhs > 1) + mexErrMsgTxt("Too many output arguments."); + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + sp = mxGetPr(prhs[1]); + + xnd = mxGetNumberOfDimensions(prhs[0]); + ynd = mxGetNumberOfDimensions(prhs[1]); + xdim = mxGetDimensions(prhs[0]); + ydim = mxGetDimensions(prhs[1]); + ndim = xnd; + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[1]); + + if(NS == 1){ + for(i=0; i<NB; i++){ + *zp++ *= *sp; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + *zp++ *= *sp++; + } + return; + } + + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + subs[i] = 0; + sx[i] = xdim[i]; + sy[i] = (i < ynd) ? ydim[i] : 1; + s[i] = sx[i] - 1; + } + + for (i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + *zp++ *= *sp; + for(i=0; i<ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/check_for_cd_arcs.m b/sourcecodes/bnt-master/BNT/potentials/check_for_cd_arcs.m new file mode 100644 index 00000000..8915d533 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/check_for_cd_arcs.m @@ -0,0 +1,23 @@ +function check_for_cd_arcs(onodes, cnodes, dag) +% CHECK_FOR_CD_ARCS Raise an error if there are any C->D links where the C node is hidden. +% check_for_cd_arcs(onodes, cnodes, dag) +% +% We cannot convert the logistic/softmax function (C->D CPD) to a Gaussian potential +% unless we use the variational approximation discussed in +% "A variational approximation for Bayesian networks with discrete and continuous latent +% variables", K. Murphy, UAI 1999. + +n = length(dag); +hnodes = mysetdiff(1:n, onodes); +chid = myintersect(cnodes, hnodes); +dnodes = mysetdiff(1:n, cnodes); +for i=chid(:)' + dcs = myintersect(children(dag, i), dnodes); + if ~isempty(dcs) + error(['hidden cts node ' num2str(i) ' has a discrete child']); + end +end + + + + diff --git a/sourcecodes/bnt-master/BNT/potentials/determine_pot_type.m b/sourcecodes/bnt-master/BNT/potentials/determine_pot_type.m new file mode 100644 index 00000000..772cc84e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/determine_pot_type.m @@ -0,0 +1,25 @@ +function pot_type = determine_pot_type(model, onodes, nodes) +% DETERMINE_POT_TYPE Determine the type of potential based on the evidence pattern. +% pot_type = determine_pot_type(model, onodes, nodes) +% +% If there are any utility nodes, pot_type = 'u' +% else +% If all hidden nodes are discrete, pot_type = 'd'. +% If all hidden nodes are continuous, pot_type = 'g' (Gaussian). +% If some hidden nodes are discrete, and some cts, pot_type = 'cg' (conditional Gaussian). +% +% nodes defaults to all nodes in graph + +nnodes = length(model.node_sizes); +if nargin < 3, nodes = 1:nnodes; end + +hnodes = mysetdiff(nodes, onodes); +if isfield(model, 'limid') %~isempty(model.utility_nodes) + pot_type = 'u'; +elseif isempty(myintersect(model.cnodes, hnodes)) + pot_type = 'd'; +elseif mysubset(hnodes, model.cnodes) + pot_type = 'g'; +else + pot_type = 'cg'; +end diff --git a/sourcecodes/bnt-master/BNT/potentials/genops.c b/sourcecodes/bnt-master/BNT/potentials/genops.c new file mode 100644 index 00000000..0416ff68 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/genops.c @@ -0,0 +1,707 @@ +/* + + GENOPS.C + Generalized arithmetic operators overloading built-in functions. + + written by Douglas M. Schwarz + schwarz@servtech.com + 26 December 1998 + Last modified: 2 April 1999 + + Copyright 1998-1999 by Douglas M. Schwarz. All rights reserved. + +*/ + + +/* + +Build MEX file by entering the appropriate command at the MATLAB prompt +(-D<name> option is equivalent to #define <name> in source file): + +mex genops.c -DPLUS_MEX -output plus +mex genops.c -DMINUS_MEX -output minus +mex genops.c -DTIMES_MEX -output times +mex genops.c -DRDIVIDE_MEX -output rdivide +mex genops.c -DLDIVIDE_MEX -output ldivide +mex genops.c -DPOWER_MEX -output power +mex genops.c -DEQ_MEX -output eq +mex genops.c -DNE_MEX -output ne +mex genops.c -DLT_MEX -output lt +mex genops.c -DGT_MEX -output gt +mex genops.c -DLE_MEX -output le +mex genops.c -DGE_MEX -output ge + +*/ + +/* This file has been formatted for a tab equal to 4 spaces. */ + +#if defined(EQ_MEX) || defined(NE_MEX) || defined(LT_MEX) || defined(GT_MEX) \ + || defined(LE_MEX) || defined(GE_MEX) +#define RELOP_MEX +#endif + +#include "mex.h" +#include "matrix.h" +#ifdef POWER_MEX +#include <math.h> +#define PI 3.141592653589793 +#endif + +bool allequal(int, const int *, const int *); +void removeZeroImag(double *, double *, int, const int *, int, mxArray **); + +#define xMat prhs[0] +#define yMat prhs[1] +#define zMat plhs[0] + +#define min(A,B) ((A) < (B) ? (A) : (B)) +#define max(A,B) ((A) > (B) ? (A) : (B)) + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + double *xrp, *xip, *yrp, *yip; +#ifndef RELOP_MEX + double *zr, *zi, *zip; +#endif + double *zrp, *zrend; + int xnd, ynd, numElements = 1; + const int *xdim, *ydim; + bool xcmplx, ycmplx; + mxClassID yclass; + int *s, ndim, *sx, *sy, i, *cpsx, *cpsy; + int *subs, *s1, *cpsx2, *cpsy2; + int ix = 0, iy = 0; + mxArray *args[3], *result[1]; +#if defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX) + double denom; +#endif +#ifdef POWER_MEX + double mag, theta, phi, magx; + int flops = 0; +#endif + + + if (nrhs != 2) + mexErrMsgTxt("Incorrect number of inputs."); + + if (nlhs > 1) + mexErrMsgTxt("Too many output arguments."); + + xnd = mxGetNumberOfDimensions(xMat); + ynd = mxGetNumberOfDimensions(yMat); + xdim = mxGetDimensions(xMat); + ydim = mxGetDimensions(yMat); + + yclass = mxGetClassID(yMat); + +/* If the built-in function in MATLAB can handle the arguments + then use that. */ + if (yclass != mxDOUBLE_CLASS || + (xnd == 2 && xdim[0] == 1 && xdim[1] == 1) || + (ynd == 2 && ydim[0] == 1 && ydim[1] == 1) || + (xnd == ynd && allequal(xnd,xdim,ydim))) + { +#ifdef PLUS_MEX + args[0] = mxCreateString("plus"); +#elif defined(MINUS_MEX) + args[0] = mxCreateString("minus"); +#elif defined(TIMES_MEX) + args[0] = mxCreateString("times"); +#elif defined(RDIVIDE_MEX) + args[0] = mxCreateString("rdivide"); +#elif defined(LDIVIDE_MEX) + args[0] = mxCreateString("ldivide"); +#elif defined(POWER_MEX) + args[0] = mxCreateString("power"); +#elif defined(EQ_MEX) + args[0] = mxCreateString("eq"); +#elif defined(NE_MEX) + args[0] = mxCreateString("ne"); +#elif defined(LT_MEX) + args[0] = mxCreateString("lt"); +#elif defined(GT_MEX) + args[0] = mxCreateString("gt"); +#elif defined(LE_MEX) + args[0] = mxCreateString("le"); +#elif defined(GE_MEX) + args[0] = mxCreateString("ge"); +#endif + args[1] = (mxArray *)xMat; + args[2] = (mxArray *)yMat; + mexCallMATLAB(1, result, 3, args, "builtin"); + mxDestroyArray(args[0]); + zMat = result[0]; + } + else /* X and Y are both N-D and different dimensionality. */ + { + ndim = max(xnd,ynd); + sx = (int *)mxMalloc(sizeof(int)*ndim); + sy = (int *)mxMalloc(sizeof(int)*ndim); + s = (int *)mxMalloc(sizeof(int)*ndim); + s1 = (int *)mxMalloc(sizeof(int)*ndim); + *(cpsx = (int *)mxMalloc(sizeof(int)*ndim)) = 1; + *(cpsy = (int *)mxMalloc(sizeof(int)*ndim)) = 1; + subs = (int *)mxMalloc(sizeof(int)*ndim); + cpsx2 = (int *)mxMalloc(sizeof(int)*ndim); + cpsy2 = (int *)mxMalloc(sizeof(int)*ndim); + for (i = 0; i < ndim; i++) + { + subs[i] = 0; + sx[i] = (i < xnd) ? xdim[i] : 1; + sy[i] = (i < ynd) ? ydim[i] : 1; + if (sx[i] == sy[i]) + s[i] = sx[i]; + else if (sx[i] == 1) + s[i] = sy[i]; + else if (sy[i] == 1) + s[i] = sx[i]; + else + { + mxFree(sx); + mxFree(sy); + mxFree(s); + mxFree(s1); + mxFree(cpsx); + mxFree(cpsy); + mxFree(subs); + mxFree(cpsx2); + mxFree(cpsy2); + mexErrMsgTxt("Array dimensions are not appropriate."); + } + s1[i] = s[i] - 1; + numElements *= s[i]; + } + + for (i = 0; i < ndim-1; i++) + { + cpsx[i+1] = cpsx[i]*sx[i]--; + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsx2[i] = cpsx[i]*sx[i]; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]); + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + xcmplx = mxIsComplex(xMat); + ycmplx = mxIsComplex(yMat); + + if (!xcmplx && !ycmplx) /* X and Y both N-D, both real. */ + { +#ifdef POWER_MEX + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX); + zrp = zr = mxGetPr(zMat); + zip = zi = mxGetPi(zMat); +#elif defined(RELOP_MEX) + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + mxSetLogical(zMat); + zrp = mxGetPr(zMat); +#else + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + zrp = mxGetPr(zMat); +#endif + xrp = mxGetPr(xMat); + yrp = mxGetPr(yMat); + zrend = zrp + numElements; + while (zrp < zrend) + { +#ifdef PLUS_MEX + *zrp++ = *xrp + *yrp; +#elif defined(MINUS_MEX) + *zrp++ = *xrp - *yrp; +#elif defined(TIMES_MEX) + *zrp++ = *xrp * *yrp; +#elif defined(RDIVIDE_MEX) + *zrp++ = *xrp / *yrp; +#elif defined(LDIVIDE_MEX) + *zrp++ = *yrp / *xrp; +#elif defined(POWER_MEX) + if (*xrp < 0.0 && *yrp != floor(*yrp)) + { + mag = pow(-*xrp,*yrp); + theta = PI * *yrp; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + *zrp++ = pow(*xrp,*yrp); + *zip++ = 0.0; + flops++; + } +#elif defined(EQ_MEX) + *zrp++ = (*xrp == *yrp); +#elif defined(NE_MEX) + *zrp++ = (*xrp != *yrp); +#elif defined(LT_MEX) + *zrp++ = (*xrp < *yrp); +#elif defined(GT_MEX) + *zrp++ = (*xrp > *yrp); +#elif defined(LE_MEX) + *zrp++ = (*xrp <= *yrp); +#elif defined(GE_MEX) + *zrp++ = (*xrp >= *yrp); +#endif + for (i = 0; i < ndim; i++) + { + if (subs[i] == s1[i]) + { + subs[i] = 0; + if (sx[i]) + xrp -= cpsx2[i]; + if (sy[i]) + yrp -= cpsy2[i]; + } + else + { + subs[i]++; + if (sx[i]) + xrp += cpsx[i]; + if (sy[i]) + yrp += cpsy[i]; + break; + } + } + } +#ifdef POWER_MEX + mexAddFlops(flops); + removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat); +#elif !defined(RELOP_MEX) + mexAddFlops(numElements); +#endif + } + else if (!xcmplx && ycmplx) /* X and Y both N-D, X real, Y complex. */ + { +#ifdef POWER_MEX + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX); + zrp = zr = mxGetPr(zMat); + zip = zi = mxGetPi(zMat); +#elif defined(RELOP_MEX) + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + mxSetLogical(zMat); + zrp = mxGetPr(zMat); +#else + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX); + zrp = mxGetPr(zMat); + zip = mxGetPi(zMat); +#endif + xrp = mxGetPr(xMat); + yrp = mxGetPr(yMat); + yip = mxGetPi(yMat); + zrend = zrp + numElements; + while (zrp < zrend) + { +#ifdef PLUS_MEX + *zrp++ = *xrp + *yrp; + *zip++ = *yip; +#elif defined(MINUS_MEX) + *zrp++ = *xrp - *yrp; + *zip++ = -*yip; +#elif defined(TIMES_MEX) + *zrp++ = *xrp * *yrp; + *zip++ = *xrp * *yip; +#elif defined(RDIVIDE_MEX) + denom = *yrp * *yrp + *yip * *yip; + *zrp++ = (*xrp * *yrp)/denom; + *zip++ = (-*xrp * *yip)/denom; +#elif defined(LDIVIDE_MEX) + *zrp++ = *yrp / *xrp; + *zip++ = *yip / *xrp; +#elif defined(POWER_MEX) + if (*yip == 0.0) + { + if (*xrp < 0.0 && *yrp != floor(*yrp)) + { + mag = pow(-*xrp,*yrp); + theta = PI * *yrp; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + *zrp++ = pow(*xrp,*yrp); + *zip++ = 0.0; + flops++; + } + } + else + { + if (*xrp < 0.0) + { + mag = pow(-*xrp,*yrp)*exp(-PI * *yip); + theta = *yip * log(-*xrp) + PI * *yrp; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + mag = pow(*xrp,*yrp); + theta = *yip * log(*xrp); + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 13; + } + } +#elif defined(EQ_MEX) + *zrp++ = (*xrp == *yrp) && (*yip == 0.0); +#elif defined(NE_MEX) + *zrp++ = (*xrp != *yrp) || (*yip != 0.0); +#elif defined(LT_MEX) + *zrp++ = (*xrp < *yrp); +#elif defined(GT_MEX) + *zrp++ = (*xrp > *yrp); +#elif defined(LE_MEX) + *zrp++ = (*xrp <= *yrp); +#elif defined(GE_MEX) + *zrp++ = (*xrp >= *yrp); +#endif + for (i = 0; i < ndim; i++) + { + if (subs[i] == s1[i]) + { + subs[i] = 0; + if (sx[i]) + xrp -= cpsx2[i]; + if (sy[i]) + { + yrp -= cpsy2[i]; + yip -= cpsy2[i]; + } + } + else + { + subs[i]++; + if (sx[i]) + xrp += cpsx[i]; + if (sy[i]) + { + yrp += cpsy[i]; + yip += cpsy[i]; + } + break; + } + } + } +#if defined(PLUS_MEX) || defined(MINUS_MEX) + mexAddFlops(2*numElements); +#elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX) + mexAddFlops(6*numElements); +#elif defined(POWER_MEX) + mexAddFlops(flops); + removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat); +#endif + } + else if (xcmplx && !ycmplx) /* X and Y both N-D, X complex, Y real. */ + { +#ifdef POWER_MEX + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX); + zrp = zr = mxGetPr(zMat); + zip = zi = mxGetPi(zMat); +#elif defined(RELOP_MEX) + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + mxSetLogical(zMat); + zrp = mxGetPr(zMat); +#else + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX); + zrp = mxGetPr(zMat); + zip = mxGetPi(zMat); +#endif + xrp = mxGetPr(xMat); + xip = mxGetPi(xMat); + yrp = mxGetPr(yMat); + zrend = zrp + numElements; + while (zrp < zrend) + { +#ifdef PLUS_MEX + *zrp++ = *xrp + *yrp; + *zip++ = *xip; +#elif defined(MINUS_MEX) + *zrp++ = *xrp - *yrp; + *zip++ = *xip; +#elif defined(TIMES_MEX) + *zrp++ = *xrp * *yrp; + *zip++ = *xip * *yrp; +#elif defined(RDIVIDE_MEX) + *zrp++ = *xrp / *yrp; + *zip++ = *xip / *yrp; +#elif defined(LDIVIDE_MEX) + denom = *xrp * *xrp + *xip * *xip; + *zrp++ = (*xrp * *yrp)/denom; + *zip++ = (-*xip * *yrp)/denom; +#elif defined(POWER_MEX) + if (*xip == 0.0) + { + if (*xrp < 0.0 && *yrp != floor(*yrp)) + { + mag = pow(-*xrp,*yrp); + theta = PI * *yrp; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + *zrp++ = pow(*xrp,*yrp); + *zip++ = 0.0; + flops++; + } + } + else + { + mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp); + theta = *yrp*atan2(*xip,*xrp); + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } +#elif defined(EQ_MEX) + *zrp++ = (*xrp == *yrp) && (*xip == 0.0); +#elif defined(NE_MEX) + *zrp++ = (*xrp != *yrp) || (*xip != 0.0); +#elif defined(LT_MEX) + *zrp++ = (*xrp < *yrp); +#elif defined(GT_MEX) + *zrp++ = (*xrp > *yrp); +#elif defined(LE_MEX) + *zrp++ = (*xrp <= *yrp); +#elif defined(GE_MEX) + *zrp++ = (*xrp >= *yrp); +#endif + for (i = 0; i < ndim; i++) + { + if (subs[i] == s1[i]) + { + subs[i] = 0; + if (sx[i]) + { + xrp -= cpsx2[i]; + xip -= cpsx2[i]; + } + if (sy[i]) + yrp -= cpsy2[i]; + } + else + { + subs[i]++; + if (sx[i]) + { + xrp += cpsx[i]; + xip += cpsx[i]; + } + if (sy[i]) + yrp += cpsy[i]; + break; + } + } + } +#if defined(PLUS_MEX) || defined(MINUS_MEX) + mexAddFlops(2*numElements); +#elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX) + mexAddFlops(6*numElements); +#elif defined(POWER_MEX) + mexAddFlops(flops); + removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat); +#endif + } + else if (xcmplx && ycmplx) /* X and Y both N-D, both complex. */ + { +#if defined(RELOP_MEX) + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + mxSetLogical(zMat); + zrp = mxGetPr(zMat); +#else + zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX); + zrp = zr = mxGetPr(zMat); + zip = zi = mxGetPi(zMat); +#endif + xrp = mxGetPr(xMat); + xip = mxGetPi(xMat); + yrp = mxGetPr(yMat); + yip = mxGetPi(yMat); + zrend = zrp + numElements; + while (zrp < zrend) + { +#ifdef PLUS_MEX + *zrp++ = *xrp + *yrp; + *zip++ = *xip + *yip; +#elif defined(MINUS_MEX) + *zrp++ = *xrp - *yrp; + *zip++ = *xip - *yip; +#elif defined(TIMES_MEX) + *zrp++ = *xrp * *yrp - *xip * *yip; + *zip++ = *xip * *yrp + *xrp * *yip; +#elif defined(RDIVIDE_MEX) + denom = *yrp * *yrp + *yip * *yip; + *zrp++ = (*xrp * *yrp + *xip * *yip)/denom; + *zip++ = (*xip * *yrp - *xrp * *yip)/denom; +#elif defined(LDIVIDE_MEX) + denom = *xrp * *xrp + *xip * *xip; + *zrp++ = (*xrp * *yrp + *xip * *yip)/denom; + *zip++ = (*xrp * *yip - *xip * *yrp)/denom; +#elif defined(POWER_MEX) + if (*xip == 0.0 && *yip == 0.0) + { + if (*xrp < 0.0 && *yrp != floor(*yrp)) + { + mag = pow(-*xrp,*yrp); + theta = PI * *yrp; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + *zrp++ = pow(*xrp,*yrp); + *zip++ = 0.0; + flops++; + } + } + else if (*xip == 0.0) + { + if (*xrp < 0.0) + { + mag = pow(-*xrp,*yrp)*exp(-PI * *yip); + theta = *yip * log(-*xrp) + PI * *yrp; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + mag = pow(*xrp,*yrp); + theta = *yip * log(*xrp); + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 13; + } + } + else if (*yip == 0.0) + { + mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp); + theta = *yrp * atan2(*xip,*xrp); + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } + else + { + magx = sqrt(*xrp * *xrp + *xip * *xip); + phi = atan2(*xip,*xrp); + mag = pow(magx,*yrp)*exp(-*yip * phi); + theta = *yip * log(magx) + *yrp * phi; + *zrp++ = mag*cos(theta); + *zip++ = mag*sin(theta); + flops += 18; + } +#elif defined(EQ_MEX) + *zrp++ = (*xrp == *yrp) && (*xip == *yip); +#elif defined(NE_MEX) + *zrp++ = (*xrp != *yrp) || (*xip != *yip); +#elif defined(LT_MEX) + *zrp++ = (*xrp < *yrp); +#elif defined(GT_MEX) + *zrp++ = (*xrp > *yrp); +#elif defined(LE_MEX) + *zrp++ = (*xrp <= *yrp); +#elif defined(GE_MEX) + *zrp++ = (*xrp >= *yrp); +#endif + for (i = 0; i < ndim; i++) + { + if (subs[i] == s1[i]) + { + subs[i] = 0; + if (sx[i]) + { + xrp -= cpsx2[i]; + xip -= cpsx2[i]; + } + if (sy[i]) + { + yrp -= cpsy2[i]; + yip -= cpsy2[i]; + } + } + else + { + subs[i]++; + if (sx[i]) + { + xrp += cpsx[i]; + xip += cpsx[i]; + } + if (sy[i]) + { + yrp += cpsy[i]; + yip += cpsy[i]; + } + break; + } + } + } +#if defined(PLUS_MEX) || defined(MINUS_MEX) + mexAddFlops(2*numElements); +#elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX) + mexAddFlops(6*numElements); +#elif defined(POWER_MEX) + mexAddFlops(flops); +#endif +#ifndef RELOP_MEX + removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat); +#endif + } + } +} + + +/*********************************************************** +* * +* Tests to see if the vectors xdim and ydim are equal. * +* * +***********************************************************/ +bool allequal(int ndim, const int *xdim, const int *ydim) +{ + int i; + bool result = true; + + for (i = 0; i < ndim; i++) + result = result && (xdim[i] == ydim[i]); + + return(result); +} + + +/****************************************************************************** +* * +* Tests to see if every imaginary element is identically zero and, if so, * +* creates a new array which is real and copies the real elements to it. * +* * +******************************************************************************/ +void removeZeroImag(double *zr, double *zi, int ndim, const int *s, + int numElements, mxArray *plhs[]) +{ + double *zrend, *ziend, *zip, *z1p, *z2p; + bool allImZero = true; + mxArray *temp; + + zip = zi; + ziend = zi + numElements; + while (zip < ziend) + { + allImZero = allImZero && (*zip++ == 0.0); + if (!allImZero) + return; + } + + temp = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + z1p = zr; + z2p = mxGetPr(temp); + zrend = z1p + numElements; + while (z1p < zrend) + *z2p++ = *z1p++; + mxDestroyArray(plhs[0]); + plhs[0] = temp; + return; +} diff --git a/sourcecodes/bnt-master/BNT/potentials/mk_initial_pot.m b/sourcecodes/bnt-master/BNT/potentials/mk_initial_pot.m new file mode 100644 index 00000000..0a52731c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/mk_initial_pot.m @@ -0,0 +1,29 @@ +function pot = mk_initial_pot(pot_type, dom, ns, cnodes, onodes) +% MK_INITIAL_POT A "initial" potential is one which has not had any evidence entered into it. +% pot = mk_initial_pot(pot_type, domain, node_sizes, cnodes, onodes) +% +% pot_type is one of 'd', 'g', 'cg' or 'u' +% domain is the set of nodes to be included in the potential. +% node_sizes(i) is the size of node i. + +switch pot_type + case 'd', + ns(onodes) = 1; + pot = dpot(dom, ns(dom)); + case 'u', + ns(onodes) = 1; + pot = upot(dom, ns(dom)); + case 'g', + ns(onodes) = 0; + pot = cpot(dom, ns(dom)); + case 'cg', + dnodes = mysetdiff(1:length(ns), cnodes); + ddom = myintersect(dnodes, dom); + cdom = myintersect(cnodes, dom); + dobs = myintersect(dnodes, onodes); + cobs = myintersect(cnodes, onodes); + ns(dobs) = 1; + ns(cobs) = 0; + pot = cgpot(ddom, cdom, ns); +end + |