diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/potentials/@scgpot')
14 files changed, 746 insertions, 0 deletions
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'); |