pa5 starter

PA3/.gitignore

@@ -1,4 +1,3 @@
 
 inf.log
 factors.fg
-*~

PA5/AssignmentToIndex.m

@@ -0,0 +1,18 @@
+% AssignmentToIndex Convert assignment to index.
+%
+%   I = AssignmentToIndex(A, D) converts an assignment, A, over variables
+%   with cardinality D to an index into the .val vector for a factor. 
+%   If A is a matrix then the function converts each row of A to an index.
+%
+%   See also IndexToAssignment.m
+
+function I = AssignmentToIndex(A, D)
+
+D = D(:)'; % ensure that D is a row vector
+if (any(size(A) == 1)),
+    I = cumprod([1, D(1:end - 1)]) * (A(:) - 1) + 1;
+else
+    I = sum(repmat(cumprod([1, D(1:end - 1)]), size(A, 1), 1) .* (A - 1), 2) + 1;
+end;
+
+end

PA5/BlockLogDistribution.m

@@ -0,0 +1,62 @@
+%BLOCKLOGDISTRIBUTION
+%
+%   LogBS = BlockLogDistribution(V, G, F, A) returns the log of a
+%   block-sampling array (which contains the log-unnormalized-probabilities of
+%   selecting each label for the block), given variables V to block-sample in
+%   network G with factors F and current assignment A.  Note that the variables
+%   in V must all have the same dimensionality.
+%
+%   Input arguments:
+%   V -- an array of variable names
+%   G -- the graph with the following fields:
+%     .names - a cell array where names{i} = name of variable i in the graph 
+%     .card - an array where card(i) is the cardinality of variable i
+%     .edges - a matrix such that edges(i,j) shows if variables i and j 
+%              have an edge between them (1 if so, 0 otherwise)
+%     .var2factors - a cell array where var2factors{i} gives an array where the
+%              entries are the indices of the factors including variable i
+%   F -- a struct array of factors.  A factor has the following fields:
+%       F(i).var - names of the variables in factor i
+%       F(i).card - cardinalities of the variables in factor i
+%       F(i).val - a vectorized version of the CPD for factor i (raw probability)
+%   A -- an array with 1 entry for each variable in G s.t. A(i) is the current
+%       assignment to variable i in G.
+%
+%   Each entry in LogBS is the log-probability that that value is selected.
+%   LogBS is the P(V | X_{-v} = A_{-v}, all X_i in V have the same value), where
+%   X_{-v} is the set of variables not in V and A_{-v} is the corresponding
+%   assignment to these variables consistent with A.  In the case that |V| = 1,
+%   this reduces to Gibbs Sampling.  NOTE that exp(LogBS) is not normalized to
+%   sum to one at the end of this function (nor do you need to worry about that
+%   in this function).
+%
+
+function LogBS = BlockLogDistribution(V, G, F, A)
+if length(unique(G.card(V))) ~= 1
+    disp('WARNING: trying to block sample invalid variable set');
+    return;
+end
+
+% d is the dimensionality of all the variables we are extracting
+d = G.card(V(1));
+
+LogBS = zeros(1, d);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% YOUR CODE HERE
+% Compute LogBS by multiplying (adding in log-space) in the correct values from
+% each factor that includes some variable in V.  
+%
+% NOTE: As this is called in the innermost loop of both Gibbs and Metropolis-
+% Hastings, you should make this fast.  You may want to make use of
+% G.var2factors, repmat,unique, and GetValueOfAssignment.
+%
+% Also you should have only ONE for-loop, as for-loops are VERY slow in matlab
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Re-normalize to prevent underflow when you move back to probability space
+LogBS = LogBS - min(LogBS);
+
+
+

PA5/CheckConvergence.m

@@ -0,0 +1,31 @@
+% CHECKCONVERGENCE Ascertain whether the messages indicate that we have converged
+%   converged = CHECKCONVERGENCE(MNEW,MOLD) compares lists of messages MNEW
+%   and MOLD.  If the values listed in any message differs by more than the 
+%   value 'thresh' then we determine that convergence has not occured and 
+%   return converged=0, otherwise we have converged and converged=1
+%
+%   The 'message' data structure is an array of structs with the following 3 fields:
+%     -.var:  the variables covered in this message
+%     -.card: the cardinalities of those variables
+%     -.val:  the value of the message w.r.t. the message's variables
+%
+%   MNEW and MOLD are the message where M(i,j).val gives the values associated
+%   with the message from cluster i to cluster j.
+%
+%
+
+function converged = CheckConvergence(mNew, mOld);
+
+thresh = 1.0e-6;
+%converged should be 1 if converged, 0 otherwise.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% YOUR CODE HERE
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+return;

PA5/ClusterGraphCalibrate.m

@@ -0,0 +1,114 @@
+% CLUSTERGRAPHCALIBRATE Loopy belief propagation for cluster graph calibration.
+%   P = CLUSTERGRAPHCALIBRATE(P, useSmart) calibrates a given cluster graph, G,
+%   and set of of factors, F. The function returns the final potentials for
+%   each cluster. 
+%   The cluster graph data structure has the following fields:
+%   - .clusterList: a list of the cluster beliefs in this graph. These entries
+%                   have the following subfields:
+%     - .var:  indices of variables in the specified cluster
+%     - .card: cardinality of variables in the specified cluster
+%     - .val:  the cluster's beliefs about these variables
+%   - .edges: A cluster adjacency matrix where edges(i,j)=1 implies clusters i
+%             and j share an edge.
+%  
+%   UseSmart is an indicator variable that tells us whether to use the Naive or Smart
+%   implementation of GetNextClusters for our message ordering
+%
+%   See also FACTORPRODUCT, FACTORMARGINALIZATION
+
+% CS228 Probabilistic Models in AI (Winter 2012)
+% Copyright (C) 2012, Stanford University
+
+function [P MESSAGES] = ClusterGraphCalibrate(P,useSmartMP)
+
+if(~exist('useSmartMP','var'))
+  useSmartMP = 0;
+end
+
+N = length(P.clusterList);
+
+MESSAGES = repmat(struct('var', [], 'card', [], 'val', []), N, N);
+[edgeFromIndx, edgeToIndx] = find(P.edges ~= 0);
+
+for m = 1:length(edgeFromIndx),
+    i = edgeFromIndx(m);
+    j = edgeToIndx(m);
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % YOUR CODE HERE
+    %
+    %
+    %
+    % Set the initial message values
+    % MESSAGES(i,j) should be set to the initial value for the
+    % message from cluster i to cluster j
+    %
+    % The matlab/octave functions 'intersect' and 'find' may
+    % be useful here (for making your code faster)
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+end;
+
+
+
+% perform loopy belief propagation
+tic;
+iteration = 0;
+
+lastMESSAGES = MESSAGES;
+
+while (1),
+    iteration = iteration + 1;
+    [i, j] = GetNextClusters(P, MESSAGES,lastMESSAGES, iteration, useSmartMP); 
+    prevMessage = MESSAGES(i,j);
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % YOUR CODE HERE
+    % We have already selected a message to pass, \delta_ij.
+    % Compute the message from clique i to clique j and put it
+    % in MESSAGES(i,j)
+    % Finally, normalize the message to prevent overflow
+    %
+    % The function 'setdiff' may be useful to help you
+    % obtain some speedup in this function
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    if(useSmartMP==1)
+      lastMESSAGES(i,j)=prevMessage;
+    end
+
+    % Check for convergence every m iterations
+    if mod(iteration, length(edgeFromIndx)) == 0
+        if (CheckConvergence(MESSAGES, lastMESSAGES))
+            break;
+        end
+        disp(['LBP Messages Passed: ', int2str(iteration), '...']);
+        if(useSmartMP~=1)
+          lastMESSAGES=MESSAGES;
+        end
+    end
+
+end;
+toc;
+disp(['Total number of messages passed: ', num2str(iteration)]);
+
+
+% Compute final potentials and place them in P
+for m = 1:length(edgeFromIndx),
+    j = edgeFromIndx(m);
+    i = edgeToIndx(m);
+    P.clusterList(i) = FactorProduct(P.clusterList(i), MESSAGES(j, i));
+end
+
+
+% Get the max difference between the marginal entries of 2 messages -------
+function delta = MessageDelta(Mes1, Mes2)
+delta = max(abs(Mes1.val - Mes2.val));
+return;
+
+

PA5/ComputeApproxMarginalsBP.m

@@ -0,0 +1,70 @@
+% COMPUTEAPPROXMARGINALSBP Computation of approximate marginals using Loopy BP
+%   M = COMPUTEAPPROXMARGINALSBP(F,E ) returns the approximate marginals over
+%       each variable in F given the evidence E.
+%   . 
+%   The Factor list F has the following fields:
+%     - .var:  indices of variables in the specified cluster
+%     - .card: cardinality of variables in the specified cluster
+%     - .val:  the cluster's beliefs about these variables
+%   - .edges: Contains indices of the clusters that have edges between them.
+%
+%   The Evidence E is a vector of length equal to the number of variables in the
+%   factors where 0 stands for unobserved and other values are an observed 
+%   assignment to that variable. It can be left empty (E=[]) if there is no evidence
+%  
+%   M should be an array of factors with one factor for each variable and 
+%   M(i).val should be filled in with the marginal of variable i.
+%
+
+
+% CS228 Probabilistic Models in AI (Winter 2012)
+% Copyright (C) 2012, Stanford University
+
+function M = ComputeApproxMarginalsBP(F,E)
+
+    % returning approximate marginals.
+
+    clusterGraph = CreateClusterGraph(F,E);
+
+    P = ClusterGraphCalibrate(clusterGraph);
+
+    N = unique([P.clusterList(:).var]);
+
+    % compute marginals on each variable
+    M = repmat(struct('var', 0, 'card', 0, 'val', []), length(N), 1);
+
+    % Populate M so that M(i) contains the marginal probability over
+    % variable i
+    for i = 1:length(N),
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % YOUR CODE HERE
+        %
+        % Populate M(i) such that M(i) contains a factor representation of
+        % the marginal proability over the variable with index i.
+        % (ie. M(i).val is the actual marginal)
+        %
+        % You may want to use the helper function 'FindPotentialWithVariable'
+        % which is defined at the bottom of this file.  Read through it
+        % to make sure you understand its functionality.
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    end
+
+
+
+return;
+
+% Helper function:
+function indx = FindPotentialWithVariable(P, V)
+
+indx = 0;
+for i = 1:length(P.clusterList),
+    if (any(P.clusterList(i).var == V)),
+        indx = i;
+        return;
+    end;
+end;
+
+return;

PA5/ComputeInitialPotentials.m

@@ -0,0 +1,79 @@
+%COMPUTEINITIALPOTENTIALS Sets up the cliques in the clique tree that is
+%passed in as a parameter.
+
+%   P = COMPUTEINITIALPOTENTIALS(C) Takes the clique tree C which is a
+%   struct with three fields:
+%   - nodes: represents the cliques in the tree.
+%   - edges: represents the adjacency matrix of the tree.
+%   - factorList: represents the list of factors that were used to build
+%   the tree. 
+%   
+%   It returns a compact form of a clique tree P that we will use through 
+%   the rest of the assigment. P is struct with two fields:
+%   - cliqueList: represents a list of cliques with appropriate factors 
+%   from factorList assigned to each clique.
+%   - edges: represents the adjacency matrix of the tree. 
+
+
+% CS228 Probabilistic Models in AI (Winter 2012)
+% Copyright (C) 2012, Stanford University
+
+function P = ComputeInitialPotentials(C)
+
+P.cliqueList = C.factorList;
+
+% number of cliques
+N = length(C.nodes);
+
+% compute assignment of factors to cliques
+alpha = zeros(length(C.factorList),1);
+
+% Setting up the cardinality 
+V = unique([C.factorList(:).var]);
+
+% Setting up the cardinality for the variables since we only get a list 
+% of factors.
+C.card = zeros(1, length(V));
+for i = 1 : length(V),
+
+    for j = 1 : length(C.factorList)
+        if (~isempty(find(C.factorList(j).var == i)))
+            C.card(i) = C.factorList(j).card(find(C.factorList(j).var == i));
+            break;
+        end
+    end
+end
+
+for i = 1:length(C.factorList),
+    for j = 1:N,
+
+        % does clique contain all variables in factor
+        if (isempty(setdiff(C.factorList(i).var, C.nodes{j}))),
+            alpha(i) = j;
+            break;
+        end;
+    end;
+end;
+
+if (any(alpha == 0)),
+    warning('Clique Tree does not have family preserving property');
+end;
+
+P.edges = C.edges;
+
+% initialize cluster potentials 
+P.cliqueList = repmat(struct('var', [], 'card', [], 'val', []), N, 1);
+
+for i = 1:N,
+    P.cliqueList(i).var = C.nodes{i};
+    P.cliqueList(i).card = C.card(P.cliqueList(i).var);
+    P.cliqueList(i).val = ones(1,prod(P.cliqueList(i).card));
+end;
+
+for i = 1:length(alpha),
+    P.cliqueList(alpha(i)) = FactorProduct(P.cliqueList(alpha(i)), C.factorList(i));    
+end;
+
+
+end
+