function wfdemo(N,L,D,Z)
% WFDEMO Demo of a simple HMM-like factorgraph with a wordfactor
%
% wfdemo(N,L,D,Z) where N is the length of the chain, L is the number
% of class labels each node can take on, Z is the feature
% dimensionality, and Z is the number of lexicon words to use.
%
% Example:
%  >> wfdemo(20,1000,100)
  
  
% Create a model with L classes and D  observation features
[Un,Ub,Uw,Uw0] = model(L,D); % Local function below

% Create some observed features for a chain plus a lexicon factor
X = features(Un,Ub,Uw,Uw0,N,L,D,Z); % 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,Uw,Uw0}, X); % Library function


% Run dense synchronous "loopy" BP with a damping factor.
fprintf('Synchronous Dense: ');
tic
  [B,M] = loopybp(G);
toc

% Run sparse synchronous loopy BP without damping.
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
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 and random
% messsage passing schedule
fprintf('Asynchronous Dense: ');
tic
  [Ba,Ma] = asyncbp(G);
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);
toc

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

%%------------------------------------------------------------------------------
% Try running the exact forward/backward pass without the lexical
% factor first, to get better initial sparsity estimates.

% Unroll a graph that does not include the lexical factors by
% leaving out their discriminants (corresponding features in X are ignored)
G0 = mkfactorgraph({Un,Ub}, X, 0); % Library function

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

fprintf('Chain-First Dense: ');
tic
  
  % Run asynchronous (forward/backward) BP with no lexicon info, and
  % keep the variable to factor messages at termination in S.
  [Ba,Ma,flag,S] = asyncbp(G0,[],[],0,order);
  
  % Run synchronous dense loopy BP on the entire graph
  [B,M] = loopybp(G,[],S);
toc


fprintf('Chain-First Sparse: ');
tic
  % Run dense asynchronous (forward/backward) BP with no lexicon info,
  % and keep the variable to factor messages at termination in S.
  [Ba,Ma,flag,S] = asyncbp(G0,[],[],0,order);
  
  % Run synchronous dense loopy BP on the entire graph
  [Bs,Ms] = sploopybp(G,1e-3,[],S);
toc


end % fgdemo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% HELPER FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Un,Ub,Uw,Uw0] = model(L,D)
% MODEL Build a model with unary and binary factors
%
% [Un,Ub,Uw,Uw0] = model(L,D) where Un and Ub and Uw0 are
% LINFACTORDISCR objects and Uw is a WORDFACTORDDISCR object.  L is
% the label cardinality (i.e., L=2 for a binary classification task),
% and D is the unary feature dimensionality. We assume no
% observations on the edges, though this is certainly possible!
  
  % Create templates for the factors  with the appropriate dimensionality
  Tn = factortemplate(L);
  Tb = factortemplate([L L]);
  Tw = factortemplate();   % No cardinalities necessary for the
                           % lexicon structural factor, allowing
                           % the  template to be used on multiple
                           % words, having different numbers of
                           % characters
  Tw0 = factortemplate(2); % Lexical indicator is binary (2-state)
  
  % 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));
  
  % Create a structural discriminant for lexicon words
  Uw = wordfactordiscr(Tw,L); % Each character takes on one of L values

  % Create a discriminant for the lexical indicator with no
  % observations and a bias feature for the indicator.
  Uw0 = linfactordiscr(Tw0,0,1); 
  
  % Set the score for indicating non-lexicon/lexicon status, respectively
  Uw0 = set(Uw0,'bias',[-abs(randn()) ; abs(randn())]);
  
end % model

function X = features(Un,Ub,Uw,Uw0,N,L,D,Z)
% FEATURES Create random observed features for a chain
%
% X = features(Un,Ub,Uw,Uw0,N,D,W) where X is an cell of factorfeature
% objects, Un is the discriminant associated with node features Ub is
% the discriminant associated with "edge" features (node pairs), Uw is
% the 0/1 discriminant associated with a lexicon indicator variable,
% and Uw0 is a bias discriminant for the lexicon variable, N is the
% chain length, L is the cardinality of each character, D is the
% dimensionality of the observed features (a scalar), and Z is the
% size of the lexicon (number of words of length N) W is the lexicon
% "matrix."
  
  % Create a cell for the N observation features, the N-1 edge
  % features, the word structural feature, and the lexical bias.
  X = cell(N+(N-1)+1+1,1);
  
  % Get the templates for the factors
  Tn = get(Un,'template');   % Node features
  Tb = get(Ub,'template');   % Edge features
  Tw = get(Uw,'template');   % Lexicon structure
  Tw0 = get(Uw0,'template'); % Lexicon prior/bias
  
  % 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

  % Create a random lexicon (L possible characters, Z words N "characters" long)
  %W = 1 + randint(L,Z,N); % This RANDINT is NOT Matlab's built-in function

  % Instead of uniformly sampling characters (as in the line above using
  % RANDINT), let's sample them from the variable-only marginals that
  % arise from the combination of the node feature discriminant
  % weights and feature vectors (so that there is some amount of
  % agreement in the support between the lexicon and the local
  % information)
  
  % First, we get the un-normalized factor tables
  M = cell2mat(cellfun(@(x) get(mkfactor(Un,x),'table'), X(1:N), ...
                       'UniformOutput',0)');

  % Since they're in log space, we need to exponentiate
  M = exp(M);
  
  % Then we normalize each variable's table to produce marginals from a
  % graph having only the node factors
  M = bsxfun(@rdivide,M,sum(M,1)); 

  % Let's add just a wee bit more randomness to make them less peaky
  M = M + rand(size(M))/L;
  M = bsxfun(@rdivide,M,sum(M,1)); 

  % Now we can  sample from these marginals to generate the lexicon
  for k=1:N
    W(:,k) = randmult(M(:,k),[Z 1]);
  end

  % Introduce the structural factor containing the lexicon, useing the
  % lexion template Tw. In this case, we only have one such feature
  % because we've modeled the problem a single word, but we can have
  % as many words/structural factors as we want over non-overlapping
  % subsets of the chain. Because the word lengths may differe, each
  % lexical factor would likely use different lexicons (word-length
  % appropriate subsets of the larger total lexicon).
  % Note that the bias variable is the N+1'th variable and it must
  % be the last entry in the variable list given to the factor feature.
  X{2*N} = factorfeature(1:(N+1), Tw, W);
  
  % Add a bias for lexicon words (we've made the lexical bias variable
  % number one more than the number of nodes in the chain), using
  % the singleton template Tw0 and no observation.
  X{2*N+1} = factorfeature(N+1,Tw0, zeros(0,1));
  
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
