Daniel@0: function [margpot, comppot] = complement_pot(pot, keep)
Daniel@0: % COMPLEMENT_POT complement means decompose of a potential into its strong marginal and 
Daniel@0: % its complement corresponds exactly to the decomposition of a probability distribution 
Daniel@0: % into its marginal and conditional
Daniel@0: % [margpot, comppot] = complement_pot(pot, keep)
Daniel@0: 
Daniel@0: % keep can only include continuous head nodes and discrete nodes
Daniel@0: % margpot is the stable CG potential of keep nodes
Daniel@0: % comppot is the stable CG potential of others in corresponds exactly to 
Daniel@0: % the discomposition of a probability distribution of its marginal and conditional
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % Calculation of the marginal requires integration over      %
Daniel@0: % all variables in csumover. Thus cheadkeep contains all     %
Daniel@0: % continuous variables in the marginal potential             %
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: %keyboard;
Daniel@0: csumover = mysetdiff(pot.cheaddom, keep);
Daniel@0: cheadkeep = mysetdiff(pot.cheaddom, csumover);
Daniel@0: 
Daniel@0: nodesizes = zeros(1, max(pot.domain));
Daniel@0: nodesizes(pot.ddom) = pot.dsizes;
Daniel@0: nodesizes(pot.cheaddom) = pot.cheadsizes;
Daniel@0: nodesizes(pot.ctaildom) = pot.ctailsizes;
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % Description of the variables in the marginal domain        %
Daniel@0: % For the calculation of a strong marginal first integration %
Daniel@0: % over all continuous variables in the head takes place.     %
Daniel@0: % The calculation of the marginal over the head variables    %
Daniel@0: % might result in a smaller or empty tail                    %
Daniel@0: % If there are no head variables, and therefore no tail      %
Daniel@0: % variables, left marginalisation over discrete variables    %
Daniel@0: % may take place                                             %
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
Daniel@0: margdom      = mysetdiff(pot.domain,keep);
Daniel@0: % margddom   = pot.ddom;
Daniel@0: margcheaddom = cheadkeep;
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % Marginalisation over discrete variables is only allowed when %
Daniel@0: % the tail is empty                                            %
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: margddom = myintersect(pot.ddom,keep);               % Discrete domain of marginal
Daniel@0: margctaildom = myintersect(pot.ctaildom,keep);       % Tail domain
Daniel@0: assert(isempty(mysetdiff(pot.ddom,margddom)) | isempty(margctaildom))  
Daniel@0: 
Daniel@0: 
Daniel@0: %margctaildom = pot.ctaildom;
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % Even if marginalisation over continuous variables is only defined %
Daniel@0: % for head variables, the marginalisation over haed-variables might %
Daniel@0: % result in a smaller tail                                          %
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: margctaildom = myintersect(pot.ctaildom,keep);
Daniel@0: 
Daniel@0: margcheadsizes = nodesizes(margcheaddom);
Daniel@0: margcheadsize = sum(margcheadsizes);
Daniel@0: margctailsizes = nodesizes(margctaildom);
Daniel@0: margctailsize = sum(margctailsizes);
Daniel@0: 
Daniel@0: compdom = pot.domain;
Daniel@0: compddom = pot.ddom;
Daniel@0: compcheaddom = csumover;
Daniel@0: compctaildom = myunion(pot.ctaildom, cheadkeep);
Daniel@0: compcheadsizes = nodesizes(compcheaddom);
Daniel@0: compcheadsize = sum(compcheadsizes);
Daniel@0: compctailsizes = nodesizes(compctaildom);
Daniel@0: compctailsize = sum(compctailsizes);
Daniel@0: 
Daniel@0: dkeep = myintersect(pot.ddom, keep);
Daniel@0: %if dom is only contain discrete node
Daniel@0: if isempty(pot.cheaddom)
Daniel@0:     dsumover = mysetdiff(pot.ddom, dkeep);
Daniel@0:     
Daniel@0:     if isempty(dsumover)
Daniel@0:         margpot = pot;
Daniel@0:         comppot = scgpot([], [], [], []);
Daniel@0:         return;
Daniel@0:     end
Daniel@0:         
Daniel@0:     
Daniel@0:     I = prod(nodesizes(dkeep));
Daniel@0:     J = prod(nodesizes(dsumover));
Daniel@0:     sum_map = find_equiv_posns(dsumover, pot.ddom);
Daniel@0:     keep_map = find_equiv_posns(dkeep, pot.ddom);
Daniel@0:     iv = zeros(1, length(pot.ddom)); % index vector
Daniel@0:     p1 = zeros(I,J);
Daniel@0:     for i=1:I
Daniel@0:         keep_iv = ind2subv(nodesizes(dkeep), i);
Daniel@0:         iv(keep_map) = keep_iv;
Daniel@0:         for j=1:J
Daniel@0:             sum_iv = ind2subv(nodesizes(dsumover), j);
Daniel@0:             iv(sum_map) = sum_iv;
Daniel@0:             k = subv2ind(nodesizes(pot.ddom), iv);
Daniel@0:             potc = struct(pot.scgpotc{k}); % violate object privacy
Daniel@0:             p1(i,j) = potc.p;
Daniel@0:         end
Daniel@0:     end
Daniel@0:     p2 = sum(p1,2);
Daniel@0:     p2 = p2 + (p2==0)*eps;
Daniel@0:     
Daniel@0:     margscpot = cell(1, I);
Daniel@0:     compscpot = cell(1, I*J);
Daniel@0:     iv = zeros(1, length(pot.ddom)); % index vector
Daniel@0:     for i=1:I
Daniel@0:         margscpot{i} = scgcpot(0, 0, p2(i));
Daniel@0:         keep_iv = ind2subv(nodesizes(dkeep), i);
Daniel@0:         iv(keep_map) = keep_iv;
Daniel@0:         for j=1:J
Daniel@0:             sum_iv = ind2subv(nodesizes(dsumover), j);
Daniel@0:             iv(sum_map) = sum_iv;
Daniel@0:             k = subv2ind(nodesizes(pot.ddom), iv);
Daniel@0:             q = p1(i,j)/p2(i);
Daniel@0:             compscpot{k} = scgcpot(0, 0, q);
Daniel@0:         end
Daniel@0:     end
Daniel@0:     
Daniel@0:     margpot = scgpot(dkeep, [], [], nodesizes, margscpot);
Daniel@0:     comppot = scgpot(pot.ddom, [], [], nodesizes,compscpot);
Daniel@0:     return;
Daniel@0: end
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % head of the potential is not empty %
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: dsize = pot.dsize;
Daniel@0: compscpot = cell(1, dsize);
Daniel@0: 
Daniel@0: fmaskh = find_equiv_posns(margcheaddom, compctaildom);
Daniel@0: fmaskt = find_equiv_posns(margctaildom, compctaildom);
Daniel@0: 
Daniel@0: fh = block(fmaskh, compctailsizes);
Daniel@0: ft = block(fmaskt, compctailsizes);
Daniel@0: 
Daniel@0: 
Daniel@0: if ~isempty(margcheaddom)
Daniel@0:     for i=1:dsize
Daniel@0:         potc = struct(pot.scgpotc{i});
Daniel@0:         q = 1;
Daniel@0:         p = potc.p;
Daniel@0:         [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(potc.A, potc.B, potc.C, margcheaddom, compcheaddom, nodesizes);
Daniel@0: 
Daniel@0:         if ~isempty(margcheaddom)
Daniel@0:             margscpot{i} = scgcpot(margcheadsize, margctailsize, p, A1, B1, C11);
Daniel@0:         else
Daniel@0:             margscpot{i} = scgcpot(margcheadsize, margctailsize, p);
Daniel@0:         end 
Daniel@0:     
Daniel@0:         if ~isempty(compcheaddom)
Daniel@0:             if ~isempty(margcheaddom)
Daniel@0:                 E = A2 - C21*pinv(C11)*A1;
Daniel@0:                 tmp1 = C21*pinv(C11);
Daniel@0:                 tmp2 = B2 - C21*pinv(C11)*B1;
Daniel@0:                 F = zeros(compcheadsize, compctailsize);
Daniel@0:                 F(:, fh) = tmp1;
Daniel@0:                 F(:, ft) = tmp2;
Daniel@0:                 G = C22 - C21*pinv(C11)*C12;
Daniel@0:             else
Daniel@0:                 E = A2;
Daniel@0:                 F = B2;
Daniel@0:                 G = C22;
Daniel@0:             end
Daniel@0:             compscpot{i} = scgcpot(compcheadsize, compctailsize, q, E, F, G);
Daniel@0:         else
Daniel@0:             compscpot{i} = scgcpot(compcheadsize, 0, q);
Daniel@0:         end
Daniel@0:         if isempty(margcheaddom)
Daniel@0:             margpot = scgpot(margddom, [], [], nodesizes, margscpot);
Daniel@0:         else
Daniel@0:             margpot = scgpot(margddom, margcheaddom, margctaildom, nodesizes, margscpot);
Daniel@0:         end
Daniel@0:     end
Daniel@0: else
Daniel@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:     % Marginalisation took place over all head variables.                               %
Daniel@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: 
Daniel@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:     % Calculate the strong marginal %
Daniel@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:     margpot = marginalize_pot(pot,keep);
Daniel@0:     mPot    = struct(margpot); 
Daniel@0:     for i =1:dsize
Daniel@0:         potc = struct(pot.scgpotc{i});  
Daniel@0:         % Get the probability of the original potential % 
Daniel@0: 	q = potc.p;
Daniel@0:          
Daniel@0:         % Get the configuration defined by the index i%
Daniel@0:         config = ind2subv(pot.dsizes,i);
Daniel@0:         
Daniel@0:         % Calculate the corresponding configuration in the marginal potential
Daniel@0:         if isempty(margpot.dsizes)
Daniel@0:             % keep == []
Daniel@0: 	    indMargPot = 1;
Daniel@0:         else
Daniel@0:             equivPos   = find_equiv_posns(dkeep,pot.ddom);
Daniel@0:             indMargPot = subv2ind(margpot.dsizes,config(equivPos));
Daniel@0:         end
Daniel@0:         % Figure out the corresponding marginal potential
Daniel@0:         mPotC = struct(mPot.scgpotc{indMargPot});
Daniel@0:         p = mPotC.p;
Daniel@0:         if p == 0
Daniel@0:             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:             % The following assignment is correct as p is only zero if q is also zero %
Daniel@0:             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:             compscpot{i} = scgcpot(compcheadsize,compctailsize,0,potc.A,potc.B,potc.C);
Daniel@0:         else
Daniel@0:             compscpot{i} = scgcpot(compcheadsize,compctailsize,q/p,potc.A,potc.B,potc.C);
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % Put all components in one potential %
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: if isempty(compcheaddom)
Daniel@0:     comppot = scgpot(compddom, [], [], nodesizes,compscpot);
Daniel@0: else
Daniel@0:     comppot = scgpot(compddom, compcheaddom, compctaildom, nodesizes,compscpot);
Daniel@0: end
Daniel@0: 
Daniel@0: