% This GCRF package was written by Jerod Weinman (jerod@acm.org). You
% are free to use it for personal or academic use. If you publish
% research that relies on this software, please kindly acknowledge it.
%
% One canonical paper you may read to understand the model the code
% is built upon is:
%
% John Laffery, Andrew McCallum, and Fernando Pereira. Conditional
% Random Fields: Probabilistic Models for Segmenting and Labeling
% Sequence Data. ICML 2001.
%
% The code works with both real-valued and binary data/feature
% vectors. Note that if you have binary data, making the data
% sparse will drastically decrease training time (orders of
% magnitude, depending on the sparsity).
%
% Here are a few examples of how to train and use a GCRF.

% Contents (search for this text to jump to the relevant part):
%% Construct a GCRF object
%% Create some training data
%% Train
%% Evaluation
%% Feature Selection
%% Instance weighting
%% Partially Labeled Data

%%
%% Construct a GCRF object
%%
% Create a GCRF with three labels, 10 node observation features and 7
% edge observation features. By default, a "bias" feature is used --
% this feature has a constant value (i.e., 1) appended to every node
% and edge feature vector and multiplied by a corresponding weight for
% each class (or pair of classes for the edge features), just like
% every other feature. 

numLabels = 3;
numNodeFeat = 10;
numEdgeFeat = 7;
isIsotropic = 1; % Horizontal and vertical edges are treated the same
isSymmetric = 1; % Left/Right and Top/Bottom are the same
defEdgeFeat = 0; % Does not append node feature vectors into edge features

C0 = gcrf(numLabels, numNodeFeat, numEdgeFeat, defEdgeFeat, isIsotropic, ...
          isSymmetric );


%%
%% Create some training data
%% 


% Create 5 instances, each with MxN = 8x4 nodes
Y = cell(1,5);

for a=1:length(Y)
  K = nodenum(8,4); % Numbers for the nodes in an 8x4 grid
  Y{a} = randint(numLabels,[8*4 1])+1; % Assign random labels
end;

% Corresponding data is represented in a special struct format with
% fields:
%  node: an M*N x numNodeFeat matrix.
%  hEdge: an M*(N-1) x numEdgeFeat matrix. 
%  vEdge: an (M-1)*N x numEdgeFeat matrix
%
% The nodes are numbered in typical column-major for the grid. The
% edges, as thought of "between" node sites, are also stored in
% column-major order (going down the rows first, then across the
% columns).
for a=1:length(Y)
  
  % Fill in node observations with random data. Note the size is
  % numFeatures x numNodes.
  X(a).node = randn(8*4,numNodeFeat);
  X(a).hEdge = randn(8*(4-1),numEdgeFeat);
  X(a).vEdge = randn((8-1)*4,numEdgeFeat);
  
end;

%%
%% Train
%%

% There are several regularizing priors available for MAP-based
% training (see the train function for more details).

% Maximum-likelihood  (train all parameters)
C1 = train(C0,Y,X,'none',[],0,[],[],[],struct('Display','final'));

% Gaussian with a standard deviation on the node weights of 10 and
% edge weights of 100.
C1 = train(C0,Y,X,'gauss',[10 1],0,[],[],[],struct('Display','final'));

% Laplace with "scale" = .1 for the nodes and 1 for the edges (in
% the dual, this means the expectations match to within that
% absolute deviation but otherwise the distribution has maximum entropy)
C1 = train(C0,Y,X,'laplace',[.1 1],0,[],[],[],struct('Display','final'));

% The original train method used an all-Matlab m-file implementation
% of belief propagation for inference. train2 is much faster (useful
% for full-resolution images), using MEX for the message-product. Go
% to inference/mex and type the command "mex calcMsgProd.c" to build
% it.

tic; 
  C1 = train(C0,Y,X,'none',[],0,[],[],[],struct('Display','final'));
toc

tic; 
  C1 = train2(C0,Y,X,'none',[],0,[],[],[],struct('Display','final'));
toc

%%
%% Evaluation
%%

% There are a few simple evaluation methods provided.

% You may get the log probability of each member of the data set,
% unfortunately the interface is outdated so you'll have to
% deconstruct X

L1 = sum(logprob(C1,Y{1},X(1).node, X(1).hEdge, X(1).vEdge));

% or the entire data set
L  = logprobd(C1,Y,X);



% Accuracy fraction (over all nodes)
% which can be calculated using MAP (default)
A = accuracy(C1,Y,X);

% or MPM
A = accuracy(C1,Y,X,1);


% You may also get the marginal probabilities for the nodes and
% edges of an instance
[Bn,Be] = marginals(C1,X(1));

% If you want these in logspace (un-exponentiated), pass a flag:
Bn = marginals(C1,X(1),1);


% Or you can simply look at the most likely joint labels yourself
P = map(C1,X(1));

% Or the most likely marginal labels
P = mpm(C1,X(1));


%%
%% Feature Selection
%%

% You can do more interesting things with training, like exclude
% certain features from use 

% Get a feature selection template matrix for the classifier
% (by default all features are excluded
S = featureSelection(C0);

% It is struct containing four matrices, one for the node weights
% and the other for the edge weights.
S

% Include all features
S.node(:) = 1;
S.hEdge(:) = 1;
S.vEdge(:) = 1;

% Completely exclude node features 2 and 5 for all 3 classes
S.node([2 5],:) = 0

% Exclude edge feature 5 for label pair (1,2)
% Because the field is symmetric, we must exclude both orderings
% (1,2) and (2,1). Because it is isotropic, we must do it for both
% h and v edges.
S.hEdge(1,2,5) = 0;
S.hEdge(2,1,5) = 0;

S.hEdge(1,2,5) = 0;
S.hEdge(2,1,5) = 0;

% Train using that feature selection
C1 = train(C0,Y,X,'gauss',[10 1],0,S);

% Look at the weights, they have zeros in the appropriate rows/features.
W = get(C1,'weights');

W.node
W.hEdge
W.vEdge
                
% You can also fix weights at non-zero values.
% Using the previous feature selection, we'll clamp those weights
% at their ML estimate:
C1 = train(C0,Y,X,'none');

W = get(C1,'weights');
W.node([1 4],:) % these will change
W.node([2 5],:) % these should stay the same

C1 = train(C1,Y,X,'gauss',10,0,S,1);

W = get(C1,'weights');
W.node([1 4],:) % these should have changed
W.node([2 5],:) % these should have stayed the same

%%
%% Instance weighting
%%

% If you want to re-weight the contribution of each training instance
% in the likelihood function (effectively modifying the so-called
% "empirical distribution") you can use a weight vector the same size
% as the labels.

% Unfortunately, I don't have this set to function at a per-node
% level, only at the whole-instance level.

% Count the third instance twice as much as the others
I = ones(size(Y));
I(3) = 2;

% Note that empty matrices can be passed for any non-scalar data to
% use defaults, and 0 can be passed for almost any scalar data.
C1 = train(C0,Y,X,'gauss',[10 1],0,[],0,I);



%%
%% Partially Labeled Data
%%

% We'll use the previously constructed data. Nodes without labels
% have ground-truth values of NaN. We'll make one node in each
% graph unlabeled.
for a=1:length(Y)
  Y{a}(floor(length(Y{a})*rand)+1) = NaN;
end;

% We can then use the function POTRAIN (partially observed
% training) with the same syntax as TRAIN. Instead of maximizing
% likelihood, it maximizes the marginal likelihood (where unlabeled
% data is marginalized). This is not convex.
C1 = potrain2(C0,Y,X,'gauss',[10 1],0,S,0,I);
