function README(N,L,D)
% FGDEMO Demo of a simple HMM-like factorgraph
%
% fgdemo(N,L,D) where N is the length of the chain, L is the number
% of class labels each node can take on, and D is the feature
% dimensionality.
%
% Example:
%  >> fgdemo(20,1000,100)
%   Synchronous Dense: Elapsed time is 53.820348 seconds.
%   Synchronous Sparse: Elapsed time is 7.821523 seconds.
%   Synchronous KL(Dense,Sparse) = 7.5556e-06
%   Asynchronous Dense: Elapsed time is 9.649072 seconds.
%   KL(Sync,Async) = 1.00533e-08
%   Asynchronous Sparse: Elapsed time is 8.076689 seconds.
%   Asynchronous KL(Dense,Sparse) = 7.47939e-06
%
% Note that while this example uses a homogeneous random field, the
% factortemplate and factordiscr classes are general and can allow
% nodes with different cardinalities and types in the same factor graph.
  
% $Id: README.m 35 2011-07-29 15:56:28Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2008, 2009, 2011
 
%    This file is part of MatLibFactorGraph.
%
%    MatLibFactorGraph is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    MatLibFactorGraph is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with MatLibFactorGraph.  If not, see <http://www.gnu.org/licenses/>.
  
% Create a model with L classes and D  observation features
[Un,Ub] = model(L,D); % Local function below
  
% Create some observed features for a chain of length 20
X = features(Un,Ub,N,D); % Local function below

% Use the factor discriminants (weights), templates, and feature
% observations to "unroll" into an actual graph. This helper factory
% method creates concrete factor objects (look-up tables associated with
% specific variables) and packages them into a factorgraph object
G = mkfactorgraph({Un,Ub}, X); % Library function


% Run dense synchronous "loopy" BP with a damping
% factor. Guaranteed to converge because we have a chain, but it
% will take longer than with a special chain-message-passing order
fprintf('Synchronous Dense: ');
tic
  [B,M] = loopybp(G);
toc

% Run sparse synchronous loopy BP Use a maximum KL divergence of 1e-3,
% corresponding to nearly exp(1-e3) = 99.9% accuracy in the messages.
fprintf('Synchronous Sparse: ');
tic
  [Bs,Ms] = sploopybp(G,1e-3);
toc

% Since M and Ms contain logs of the node marginals (which should be exact
% for M), we can compute  the KL divergence from Ms to M in bits to
% measure the error.

fprintf('Synchronous KL(Dense,Sparse) = %g\n',...
        sum(cellfun(@(me,ms) sum( exp(me).*(me-ms) ), M,  Ms)) )/log(2);

% Run dense asynchronous "loopy" BP with no damping
% factor. Specify the order of the message passing so that it
% should converge in only two iterations, a forward and backward pass

% Get the chain ordering of the factors
order = factororder(N); % Helper below

fprintf('Asynchronous Dense: ');
tic
  [Ba,Ma] = asyncbp(G,[],[],0,order);
toc

fprintf('KL(Sync,Async) = %g\n',...
        sum(cellfun(@(me,ma) sum( exp(me).*(me-ma) ), M,  Ma)) )/log(2);

% Run sparse asynchronous loopy BP. With small chains and
% statespaces, this will take longer because it will have
% extraneous sorts and intermediate updates.
fprintf('Asynchronous Sparse: ');
tic
  [Bas,Mas] = spasyncbp(G,1e-3,[],[],0,order);
toc

fprintf('Asynchronous KL(Dense,Sparse) = %g\n',...
        sum(cellfun(@(ma,mas) sum( exp(ma).*(ma-mas) ), Ma,  Mas)) )/log(2);


end % fgdemo


function [Un,Ub] = model(L,D)
% MODEL Build a model with unary and binary factors
%
% [Un,Ub] = model(L,D) where Un and Ub are LINFACTORDISCR objects. L
% is the label cardinality (i.e., L=2 for a binary classification
% task), and D is the unary feature dimensionality. We assume no
% observationson the edges, though this is certainly possible!
  
  % Create templates for the factors with the appropriate dimensionality
  Tn = factortemplate(L);
  Tb = factortemplate([L L]);
  
  % Create a linear discriminant for the unary factor with
  % dimensionality D and a bias feature
  Un = linfactordiscr(Tn,D,1);
  
  % Use random weights for linear discriminant (one weight vector
  % for each class)
  Un = set(Un,'weights', randn(L,D));
  
  % Set a random bias vector (implicity multiplied with a feature
  % of value 1)
  Un = set(Un,'bias',randn(L,1));
  
  % Create a linear discriminant for the binary factor with no
  % observations and bias features for the class labels (i.e., a
  % character bigram)
  Ub = linfactordiscr(Tb,0,1);

  % Set a random bias vector (implicity multiplied with a feature
  % of value 1) 
  Ub = set(Ub,'bias',randn(L^2,1));
end % model

function X = features(Un,Ub,N,D)
% FEATURES Create random observed features for a chain
%
% X = features(In,Ub,Tn,Tb,N) where X is an cell of factorfeature
% objects, Tn is a unary factortemplate object to associate with
% the features, Tb is the binary factor template to connect the
% features with, N is the chain length, and D is the dimensionality
% of the observed features (a scalar).
  
  % Get the templates for the factors
  Tn = get(Un,'template');   % Node features
  Tb = get(Ub,'template');   % Edge features

  % Create a cell for the N observation features and the N-1 edge features
  X = cell(N+(N-1),1);
  
  % Loop to create N observed features
  for k=1:N
    % Create the unary factor feature that has a single variable
    % number (k), uses the unary template Tn, and has a random observation
    X{k} = factorfeature([k], Tn, randn(D,1));
  end % node loop
  
  
  % Loop to create N-1 edge features
  for k=1:N-1
    % Create the binary factor feature that has a neighboring pair of variables
    % (k and k+1), uses the binary template Tb, and has no observation
    X{N+k} = factorfeature([k k+1], Tb, zeros(0,1) );
  end % edge loop
  
end % features


function F = factororder(N)
% ORDER Topological order of the factor indices of a chain 
  
% Nodes were given first, then edges, so we must interleave them
  F = zeros(1,N+(N-1));
  m = 1;
  for k=1:N 
    F(m) = k; % Node factor index

    if (k ~= N) % No edge to right of last node
      F(m+1) = N+k; % Factor index of edge to right of node
    end
    
    m = m+2;
  end
  
  end %factororder
