function [C,fval,exitflag,output] = train2(C,Y,X,P,R,N,S,F,I,options)
% TRAIN Set the parameters using MAP with synchronous loopybp
%
% TRAIN(C,Y,X,[P],[R],[N],[S],[F],[I],[options]) where C is the gcrf, Y is a
% cell of node label matrices, X is an array of "grid graph feature"
% structs, P specifies `the regularizing prior
% {'gauss','laplace','none'}, R is the regularization
% parameter(s). N=1 means to train only node parameters (disregarding
% adjacencies) N=3 means to train edge parameters with node params
% fixed. N=0 trains all parameters. N=2 trains node parameters with
% edge params fixed.  S is a feature selection. F indicates whether
% unselected features have zero weights (0) or are frozen at current
% values. I are instance weights.
  
% R is [nodeParam edgeParam]


% Jerod Weinman
% jerod@acm.org
    
% $Id: train2.m 53 2012-12-28 15:17:20Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2004, 2008, 2009
 
%    This file is part of Matlab GridCRF.
%
%    Matlab GridCRF 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.
%
%    Matlab GridCRF 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 Matlab GridCRF.  If not, see <http://www.gnu.org/licenses/>.

% GCRF.TRAIN
  
% Handle arguments and set defaults

  % Regularization method
  if (nargin>=4)
	
	if strcmpi(P,'gauss')
	  P = 2;
	elseif strcmpi(P,'laplace');
	  P = 1;
	elseif strcmpi(P,'none');
	  P = 0;
	else
	  error(sprintf('Unrecognized regularization method ''%s''',P));
	end;
  else
	P = 2;
  end;
  
  % Regularization method's parameters (sigma if gaussian, alpha if laplacian)
  if (nargin<5) || isempty(R)
	switch P
	 case 0
	  R = [0 0];
	 case 1
	  R = [.1 .01];
	 case 2
	  R = [10 1];
	end;
  else
	if any(R<0)
	  error('Negative regularization parameters are not permitted');
	end;
	
	if length(R)==1
	  R = [R R];
	end;
	
  end;
  
  % Which features to train (default all)
  if (nargin<6)
	N = 0;
  elseif N~=0 && N~=1 && N~=2 && N~=3
	error('Unknown argument for parameter training type N');
  end
  
  % Feature selection
  if (nargin<7)
	S = [];
  elseif ~isempty(S)
	np = sumv(S.node) + sumv(S.hEdge) + sumv(S.vEdge);

	if (C.defEdge)
	  np = np + sumv(S.dhEdge) + sumv(S.dvEdge);
	end;
	
	disp(['numParameters = ' num2str(np)]);
	
  end;
  
  % Whether to fix unselected parameters or zero them
  if nargin<8
	F = 0;
  end;
  
  % Instance weights
  if nargin<9 | isempty(I)
	I = ones(1,length(X));
  elseif length(I)~=length(X)
	error('Instance weights length do not match number of instances');
  elseif any(I<0)
	error('Instance weights must be positive');
  end;
  
  % Optimizer parameters/options
  if nargin<10
	options = [];
  end;
	
  
  if length(Y) ~= length(X)
	error('Labels and feature graphs lengths do not match');
  end;
  
  % Make sure persistent variables are cleared in these functions
  paramsToWeights();
  weightsToParams();
  
  
  % Pass in weights to set dimensions:
  % "weights" = struct, "params" = vector for optimizer
  paramsToWeights(C.weights, C.defEdge, C.isotropic, C.symmetric);
  
  Wv = weightsToParams(C.weights, C.defEdge, C.isotropic, C.symmetric);

  disp(sprintf(['Training gcrf:\n\tdefEdge = %d\n\tisotropic = %d\n\t' ...
				'symmetric = %d\n\tnumLabels = %d\n\tnumNodeParams = %d\n\t' ...
				'numEdgeParams = %d\n\tprior = %d\n\tparams = [%f %f]'],...
			   C.defEdge,C.isotropic,C.symmetric, C.numLabels, ...
			   size(C.weights.node,1)-1,size(C.weights.hEdge,3)-1, ...
			   P,R(1), R(2)));

  
  % Pass in training data and params to the value function
  value(C.weights,Y,X,P,R,N,S,F,I,C);
  
  
  % Get hook to the value function
  F = @value;
  
  
  
  % Get the initial value
  [V,G] = value(Wv); 
  
  if 1
	  
	exitflag = -3;
	  
	if isempty(options)
	  options = struct('Display','iter',...
					   'MaxIter',1000,...
					   'TolFun',1e-4,...
					   'TolX', 1e-10,...
					   'TolGrad',1e-4);
	end;		
	  
	while (exitflag == -3)
	  [Wv,fval,exitflag,output] = lbfgs( F, Wv, options);
		
	  fval
	  exitflag
	  output
	end;
	  
  else
	% Keep this around for the useful derivative checks
	options = optimset('Display','iter',...
					   'MaxIter',1,...
					   'GradObj','on',...
					   'Hessian','off',... %'TolX',1e-10,...
					   'LargeScale','on',...
					   'DerivativeCheck','on');
					   
					   
	[Wv,fval,exitflag,output,grad] = fminunc(F, Wv, options);
	  
	fval
	exitflag
	output

  end;

  
  % Clear the persistent vars
  weightsToParams();
	 
  V
  
  % Set the weights we just received from the optimizer
  C.weights = paramsToWeights(Wv);
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = paramsToWeights(P, varargin)
% paramsToWeights transforms a parameter vector to a weight struct
  
  persistent numNodeFeat numHorizFeat numVertFeat numLabels ...
	  ndxHorizFeat ndxVertFeat tieHorizFeat tieVertFeat D I S;

  if nargin==0
	clear numNodeFeat numHorizFeat numVertFeat numLabels ...
		ndxHorizFeat ndxVertFeat tieHorizFeat tieVertFeat D I S;
	return;
  end;
  
  % If passed the actual weights struct, 
  if isa(P,'struct')

	D = varargin{1}; % Has Default edge features
	I = varargin{2}; % Is Isotropic
	S = varargin{3}; % Is Symmetric
	
	
	% initialize persistent vars
	numNodeFeat = size( P.node, 1);
	numHorizFeat = size( P.hEdge, 3);
	numVertFeat = size( P.vEdge, 3);
	
	numLabels = size( P.node, 2);
	

	if S % symmetric
	  
	  % since edge features are symmetric in the label, we take the
	  % upper-triangular (label1[row] < label2[col]) weights and
	  % stick them in the lower-triangular parts
	
	  LL =  triu(ones(numLabels));
	  
	  ndxHorizFeat = find(LL(:,:,ones(1,numHorizFeat)));
	  
	  
	  [c1,c2,ff] = ind2sub(size(P.hEdge),ndxHorizFeat);
	  
	  tieHorizFeat = sub2ind(size(P.hEdge), c2, c1, ff);
	  
	  if ~I
		ndxVertFeat = find(LL(:,:,ones(1,numVertFeat)));
		
		[c1,c2,ff] = ind2sub(size(P.vEdge),ndxVertFeat);
		
		tieVertFeat = sub2ind(size(P.vEdge), c2, c1, ff);
	  end;
	end;
	
	% and return -- we're only setting persistent vars, now we're done
	return;
  end;
  

  nnw = numNodeFeat * numLabels; % Number of node weights

  
  % Pull the node weights out of the parameter vector, and reshape appropriately
  W.node = reshape( P(1:nnw), [numNodeFeat numLabels]);
  

  
  if S % symmetric
	
	ndx = nnw+(1:length(ndxHorizFeat));
	
	W.hEdge = zeros([ numLabels numLabels numHorizFeat]);
	
	W.hEdge(ndxHorizFeat) = P(ndx);
	W.hEdge(tieHorizFeat) = P(ndx);
	
	if I % symmetric and isotropic (take symmetric horiz features)
	  
	  W.vEdge = zeros([ numLabels numLabels numVertFeat]);

	  W.vEdge(ndxHorizFeat) = P(ndx);
	  W.vEdge(tieHorizFeat) = P(ndx);
	else % symmetric and anisotropic
	  
	  W.vEdge = zeros([ numLabels numLabels numVertFeat]);
	  
	  ndx = ndx(end) + (1:length(ndxVertFeat));

	  W.vEdge(ndxVertFeat) = P(ndx);
	  W.vEdge(tieVertFeat) = P(ndx);
	end;
	
  else % asymmetric
	
	ndx = nnw+(1:(numLabels*numLabels*numHorizFeat));
	
	W.hEdge = reshape( P(ndx), [ numLabels numLabels numHorizFeat]);
  
	if I % asymmetric and isotropic -- this does't really make sense

	  W.vEdge = W.hEdge;
	
	else % asymmetric and anisotropic -- the "common case"?
	  
	  ndx = ndx(end) + (1:(numLabels*numLabels*numHorizFeat));

	  W.vEdge = reshape( P(ndx), [ numLabels numLabels numHorizFeat]);
	
	end;

  end;
  
  if D % default edge features used
	
	if S % symmetric
	  ndx = ndx(end) + (1:(numLabels*numLabels*(numNodeFeat-1)));
	  
	  W.dhEdge = reshape( P(ndx), [numLabels numLabels numNodeFeat-1]);
	  W.dhEdge(:,:,:,2) = permute(W.dhEdge(:,:,:,1), [2 1 3 4]);
	  
	  if I % symmetric and isotropic (take symmetric horiz features)
		W.dvEdge = W.dhEdge;
	  else % symmetric and anisotropic 
		
		ndx = ndx(end) + (1:(numLabels*numLabels*(numNodeFeat-1)));

		W.dvEdge = reshape( P(ndx), [numLabels numLabels numNodeFeat-1]);
		W.dvEdge(:,:,:,2) = permute(W.dvEdge(:,:,:,1), [2 1 3 4]);
	  end;
	  
	else % asymmetric
	  
	  ndx = ndx(end) + (1:(numLabels*numLabels*(numNodeFeat-1)*2));
	  
	  W.dhEdge = reshape( P(ndx), [numLabels numLabels numNodeFeat-1 2]);
	  
	  if I % asymmetric and isotropic
		W.dvEdge = W.dhEdge;
	  else % asymmetric and anisotropic
		ndx = ndx(end) + (1:(numLabels*numLabels*(numNodeFeat-1)*2));
		
		W.dvEdge = reshape( P(ndx), [numLabels numLabels numNodeFeat-1 2]);
	  end;
	  
	end;
  else
	W.dhEdge = zeros([numLabels numLabels numNodeFeat-1 2]);
	W.dvEdge = zeros([numLabels numLabels numNodeFeat-1 2]);
  end;
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [P,N] = weightsToParams(W,varargin)
% weightsToParams transforms a weight struct to a parameter vector
  
  persistent D I S ndxHorizFeat ndxVertFeat
  

  if nargin==0
	clear D I S ndxHorizFeat ndxVertFeat;
	return;
  elseif nargin>1
	D = varargin{1};
	I = varargin{2};
	S = varargin{3};
  end;
  
  if S & isempty(ndxHorizFeat)
	
	nL = size(W.node,2);
	nH = size(W.hEdge,3);
	nV = size(W.vEdge,3);
	
	LL =  triu(ones(nL));
	  
	ndxHorizFeat = find(LL(:,:,ones(1,nH)));

	if ~I
	  ndxVertFeat = find(LL(:,:,ones(1,nV)));
	end;
  end;

  
  P = W.node(:)';
  
  N = length(P);
  
  if S
	P = [P W.hEdge(ndxHorizFeat)'];
	
	if ~I
	  P = [P W.vEdge(ndxVertFeat)'];
	end;
  else
	P = [P W.hEdge(:)'];
	
	if ~I
	  P = [P W.vEdge(:)'];
	end;
  end;
  
  if D
	if S
	  P = [P mkrvector(W.dhEdge(:,:,:,1))];
	  
	  if ~I
		P = [P mkrvector(W.dvEdge(:,:,:,1))];
	  end;
	else
	  P = [P W.dhEdge(:)'];
	  
	  if ~I
		P = [P W.dvEdge(:)'];
	  end;
	end;
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [V,G] = value(varargin)
  
% A few constant values for clarity
  ALL = 0;
  NODE_ONLY = 1;
  EDGE_FIXED = 2;
  NODE_FIXED = 3;
  
  persistent Y X P R N S F I C Gobs Wold Vold Gold ;
  
  switch length(varargin) 
   case 1
	W = varargin{1}; % weights vector
   case 10;
	
	% Set the values of the persistent variables
	W = varargin{1}; % weights STRUCT
	Y = varargin{2}; % labels
	X = varargin{3}; % data
	P = varargin{4}; % regularization prior
	R = varargin{5}; % regularization prior params
	N = varargin{6}; % which parameters
	S = varargin{7}; % feature selection
	F = varargin{8}; % freeze unselected weights
	I = varargin{9}; % instance weights
	C = varargin{10};% the gcrf

	% Store observations
	Gobs = calcObsConstraints(W,Y,X,N,I,C.defEdge,C.isotropic,C.symmetric);
		
	Wold = [];
	% And exit
	return;
	
   otherwise
	error('Incorrect arguments');
  end;
    
  % Use cached value and gradient
  if ~isempty(Wold) & max(abs(Wold-W))<5*eps
	V = Vold;
	G = Gold;
	disp(sprintf('cached value = %g\tmax(|w-w''|) = %g',V, max(abs(Wold-W))));
	return;
  end;

  %if nargout==2
	doGrad = 1; % if these conditionals are uncommented, caching
                % will break
  %else
  %	doGrad = 0;
  %end;



  Ws = paramsToWeights(W); % Get the weight struct

  disp(sprintf('|weights| = %g\t|node| = %g\t|edge| = %g',...
			   norm(W,1),norm(Ws.node(:),1),norm(W,1)-norm(Ws.node(:),1)));
  
  % Feature selection
  if ~isempty(S) & ~F
	for f=fieldnames(S)'
	  Ws.(f{:}) = Ws.(f{:}) .* S.(f{:});
	end;
  end;
  
  V = 0;

  if doGrad
	% Gradient = Observation - Expectation
	G = Gobs;
  end;
  
  for a=1:length(Y)

	M = unrollgridflat(Ws,X(a).node, X(a).hEdge,X(a).vEdge, C.defEdge);

    K = nodenum(gridsize(X(a).node,X(a).hEdge));
    
	%tic;
    if N==NODE_ONLY
      B = gridbeliefs(K,M,[],0,1);
    else
      [B,flag] = gridsyncbp(K, M, 0,[],1);
      %[B,flag] = gridbp(K, M, 0,[],1);
      %disp(sprintf('BP time: %f',toc));
	end;

	%if (flag>0) % did BP converge?
	  vv = gridlogprob(K,Y{a}', B ,1);
	%else

	  % BP did not converge, so these are bad or useless params. Give an
      % infinite value so the optimizer will back off.
	  
	  %vv = -Inf;
	%end;
	
	% BP returns beliefs in logspace. exponentiate them so we can
    % use them to weight the features for the expectations.
    B.node = exp(B.node);
    B.hEdge = exp(B.hEdge);
    B.vEdge = exp(B.vEdge);
    
	% Add (weighted) instance value/logprob to total
	V = V + I(a)*vv;
	
	
	if isinf(vv)
	  warning(sprintf('Infinite value for instance %d',a));
	  
	  % Return infinite value. 
	  V = Inf;

	  % Gradient will be meaningless anyway, so use regularization grad
	  [rV,rG] = regValue(Ws,P,R,N);
	  G = - weightsToParams(rG);
	  return; 
	elseif isnan(vv)
	  warning(sprintf('NaN value for instance %d',a));

	  % Return NaN value. 
	  V = NaN;
	  % Gradient will be meaningless anyway, so use regularization grad
	  [rV,rG] = regValue(Ws,P,R,N);
	  G = - weightsToParams(rG);

	  return;
	elseif 0 && (vv>eps)
	  warning(sprintf('Positive log likelihood (instance %d) value %f', a,vv));
	  
	  % I don't know what this means. It COULD be an error in
      % asyncbp.m or logprob.m, or it could result from the fact we
      % have an approximation.
	  
	  % Even if it makes sense (not sure yet) let's not allow
      % positive likelihoods by returning an infinite value
	 
	  V = Inf;
	  % Gradient will be meaningless anyway, so use regularization grad
	  [rV,rG] = regValue(Ws,P,R,N);
	  G = - weightsToParams(rG);
	  return;
	end;
					  
	
	if doGrad
	  
	  Gexp = calcObsExpectations(C,X(a),I(a),B,N);
	  
	  for f=fieldnames(Gexp)'
		G.(f{:}) = G.(f{:}) - Gexp.(f{:});
	  end;
	  
	end;
	
	
  end;


  disp(sprintf('raw value = %g',V));
  
  [rV,rG] = regValue(Ws,P,R,N);
  
  
  V = V + rV;

  disp(sprintf('regularized value = %g',rV));

  
  if doGrad
	% Add the prior/regularization to gradient
	if N~=NODE_FIXED
	  G.node = G.node + rG.node;
	end;
	
	if N~=EDGE_FIXED & N~=NODE_ONLY
	  G.hEdge = G.hEdge + rG.hEdge;
	  G.vEdge = G.vEdge + rG.vEdge;
	  
	  if C.defEdge
		G.dhEdge = G.dhEdge + rG.dhEdge;
		G.dvEdge = G.dvEdge + rG.dvEdge;
	  end;
	end;
	
	% Zero the gradient of unselected features so the optimizer
    % won't change those params
	if ~isempty(S)
	  for f=fieldnames(G)'
		G.(f{:}) = G.(f{:}).*S.(f{:});
	  end;
	end;

	% Flatten!
	G = weightsToParams(G); 
	
	disp(sprintf('|grad| = %g',norm(G,1)));
  end;
  
  V = -V;  % Since the optimization toolbox does minimization
  G = -G;

  % Cache
  Wold = W;
  Vold = V;
  
  if doGrad
	Gold = G;
  else
	Gold = NaN;
  end;
  
  
  disp(sprintf('total value = %f',V));
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [V,G] = regValue(W,P,R,N)
% REGVALUE Regularization value and gradient
  
  ALL = 0;
  NODE_ONLY = 1;
  EDGE_FIXED = 2;
  NODE_FIXED = 3;
  UNSL_FIXED = 4;
  
  
  % Flatten to make our lives easy, and grab the number of node
  % parameters, which will be in the front.
  
  [Wv,np] = weightsToParams(W);
  Gv = zeros(size(Wv));

  % Regularization
  switch P
   case 0 % No prior p(w) \propto 1
	V = 0; 
   case 1 % Laplacian p(w) \propto exp(-a|w|)
	
	V = 0;

	if N~=NODE_FIXED
	  V =  V - sum(abs(Wv(1:np)))*R(1);

	  Gv(1:np) =  - sign(Wv(1:np))*R(1);
	end;
	
	if N~=EDGE_FIXED & N~=NODE_ONLY
	  V = V - sum(abs(Wv(np+1:end)))*R(2);
	
	  Gv(np+1:end) =  - sign(Wv(np+1:end))*R(2);
	end;
	
   case 2 % Gaussian p(w) \propto exp(-1/2*w^2/s^2)

	V = 0;

	if N~=NODE_FIXED
	  V = - sum(Wv(1:np).^2)./(2*R(1)^2);
	
	  Gv(1:np) = - Wv(1:np)./R(1)^2;
	end;
	
	if N~=EDGE_FIXED & N~=NODE_ONLY
	  V = V - sum(Wv(np+1:end).^2)./(2*R(2)^2);

	  Gv(np+1:end) = - Wv(np+1:end)./R(2)^2;
	end;
	
   otherwise
	error('Unknown regularization');
  end;
  
  % Return to struct format
  G = paramsToWeights(Gv);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F = calcObsConstraints(W,Y,X,N,I,D,P,S);

% Calculate feature observations (like expectations, but with delta
% functionson the labels)
  
  ALL = 0;
  NODE_ONLY = 1;
  EDGE_FIXED = 2;
  NODE_FIXED = 3;
  UNSL_FIXED = 4;
  
  disp('Calculating constraints ... ');
  
  L = size(W.node,2); % Number of node labels

  F = W;
  
  for f=fieldnames(F)'
	F.(f{:})(:) = 0;
  end;
  
  % Do this once at the beginning, and once at the end, so we don't
  % have to do it many times in a loop.
  F.node = F.node';
  
  for a=1:length(Y)
	
	yy = Y{a};
	xx = X(a);
	
	sz = gridsize(xx.node,xx.hEdge); % get grid dimensions
	K = nodenum(sz); % number the nodes for each indexing of yy
	
	
	% Bn   LxN
	% Be{e} LxL
	
	if N~=NODE_FIXED
		
	  Bn = sparse(L,length(yy));

	  ndx = sub2ind(size(Bn), yy(:)', 1:size(Bn,2));
		
	  Bn(ndx) = 1;
	  
	  F.node = F.node + ...
			   I(a) * [Bn*xx.node Bn*ones(size(xx.node,1),1)];
	  
	end;
	
	if N~=NODE_ONLY & N~=EDGE_FIXED
		
	  % Horizontal Features
	  ff = zeros(size(F.hEdge));

	  % Left nodes for horizontal edges
	  [R,C] = ind2sub([sz(1) sz(2)-1],1:size(xx.hEdge,1));

	  for e=1:length(R)

		l0 = yy(K(R(e),C(e)));
		l1 = yy(K(R(e),C(e)+1));
		
		if S % symmetric
		  [l0,l1] = order(l0,l1);
		end;
		  
		ff(l0,l1,:) = ff(l0,l1,:) + permute([xx.hEdge(e,:) 1],[1 3 2]);
		
	  end;
	  
	  F.hEdge = F.hEdge + I(a) * ff;
	  
	  % Vertical Features
	  ff = zeros(size(F.vEdge));
	  
	  % Top nodes for vertical edges
	  [R,C] = ind2sub([sz(1)-1 sz(2)], 1:size(xx.vEdge,1));

	  for e=1:length(R)

		l0 = yy(K(R(e),C(e)));
		l1 = yy(K(R(e)+1,C(e)));
		
		if S % symmetric
		  [l0,l1] = order(l0,l1);
		end;		  
		
		ff(l0,l1,:) = ff(l0,l1,:) + permute([xx.vEdge(e,:) 1],[1 3 2]);
	  end;

	  
	  F.vEdge = F.vEdge + I(a) * ff;
	  
	  
	  if P % isotroPic
		% Smoosh horizontal and vertical together (into horizontal
        % slots) and call it a day (zero vertical slots)
		F.hEdge = F.hEdge + F.vEdge;
		F.vEdge(:) = 0;
	  end;
	  
	  if D % default edge features
	
		% Horizontal Features
		ff = zeros(size(F.dhEdge));

		% Left nodes for horizontal edges
		[R,C] = ind2sub([sz(1) sz(2)-1],1:size(xx.hEdge,1));
		
		for e=1:length(R)
		  
		  n0 = K(R(e),C(e));
		  n1 = K(R(e),C(e)+1);
		  
		  l0 = yy(n0);
		  l1 = yy(n1);
		  
		  ff(l0,l1,:,1) = ff(l0,l1,:,1) + permute(xx.node(n0,:),[1 3 2]);
		  ff(l0,l1,:,2) = ff(l0,l1,:,2) + permute(xx.node(n1,:),[1 3 2]);
		  
		end;

		F.dhEdge = F.dhEdge + I(a) * ff;
		
		% Vertical Features
		ff = zeros(size(F.dvEdge));

		% Top nodes for vertical edges
		[R,C] = ind2sub([sz(1)-1 sz(2)], 1:size(xx.vEdge,1));

		for e=1:length(R)
		  
		  n0 = K(R(e),C(e));
		  n1 = K(R(e)+1,C(e));

		  l0 = yy(n0);
		  l1 = yy(n1);
		
		  ff(l0,l1,:,1) = ff(l0,l1,:,1) + permute(xx.node(n0,:),[1 3 2]);
		  ff(l0,l1,:,2) = ff(l0,l1,:,2) + permute(xx.node(n1,:),[1 3 2]);
		end;
		
		
		F.dvEdge = F.dvEdge + I(a) * ff;
		
		if S % symmetric

		  % Need to reverse the label ordering on the second node's
          % features, and add it to the first node features (more smooshing)
		  F.dhEdge(:,:,:,1) = F.dhEdge(:,:,:,1) + ...
			  permute(F.dhEdge(:,:,:,2),[2 1 3 4]);
		  
		  F.dhEdge(:,:,:,2) = 0;
		  

		  F.dvEdge(:,:,:,1) = F.dvEdge(:,:,:,1) + ...
			  permute(F.dvEdge(:,:,:,2),[2 1 3 4]);
		  
		  F.dvEdge(:,:,:,2) = 0;
		  
		end;
		
		if P % isotroPic
			 % Smoosh horizontal and vertical together (into horizontal
			 % slots) and call it a day (zero vertical slots)
		  F.dhEdge = F.dhEdge + F.dvEdge;
		  F.dvEdge(:) = 0;
		end;
		
	  end; % default edge features

	end; % edge params
	
  end; % instance loop

  F.node = F.node';
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G = calcObsExpectations(C,X,I,B,N)
  
  ALL = 0;
  NODE_ONLY = 1;
  EDGE_FIXED = 2;
  NODE_FIXED = 3;
  UNSL_FIXED = 4;
  
  
  G = C.weights;
  
  for f=fieldnames(G)'
	G.(f{:})(:) = 0;
  end;
	
  % Bn   LxN
  % Be{e} LxL

  sz = gridsize(X.node,X.hEdge);
  L = C.numLabels;
  
  if (N~=NODE_FIXED) % fix node weights -> zero gradient
	
	G.node = G.node + ...
			 I * [B.node*X.node  B.node*ones(size(B.node,2),1)]';	  
	
  end;
  
  if (N~=EDGE_FIXED && N~=NODE_ONLY)
	
	K = nodenum(sz);
	
    
	% Horizontal Features
	gg = zeros(size(G.hEdge));
	fl = size(G.hEdge,3);
	
	for e=1:size(B.hEdge,1)
		
		be = B.hEdge(e,:)'; % L^2x1
		
		if C.symmetric
		  % Don't distinguish the ordering of the labels
          be = reshape(be,[L L]);
		  be = triu(be) + tril(be,-1)';
          be = be(:); % vectorize: L^2x1
		end;
		
		bex = [be*X.hEdge(e,:) be];
		
		gg = gg + reshape(bex,[L L fl]);
		
	end;
	
	G.hEdge = G.hEdge + I*gg;
	
	% Vertical Features
	gg = zeros(size(G.vEdge));
	fl = size(G.vEdge,3);

	for e=1:size(B.vEdge,1)
		
		be = B.vEdge(e,:)'; % L^2x1
		
		if C.symmetric
		  % Don't distinguish the ordering of the labels
          be = reshape(be,[L L]);
		  be = triu(be) + tril(be,-1)';
          be = be(:); % vectorize: L^2x1
		end;
		
		bex = [be*X.vEdge(e,:) be];
		
		gg = gg + reshape(bex,[L L fl]);
		
	end;
	
	G.vEdge = G.vEdge + I*gg;
	
	if C.isotropic
	  G.hEdge = G.hEdge + G.vEdge;
	  G.vEdge(:) = 0;
	end;
	
	if C.defEdge
	  
	  % Horizontal Features
	  szd = size(G.dhEdge);

	  fl = size(X.node,2);
	  
	  gg1 = zeros(szd(1:3));
	  gg2 = zeros(szd(1:3));
	  
	  [R,M] = ind2sub([sz(1) sz(2)-1],1:size(X.hEdge,1));
	
	  for e=1:length(R)
		
		k1 = K(R(e),M(e));
		k2 = K(R(e),M(e)+1);
		
		be = B.hEdge(e,:)'; % L^2x1
		
		gg1 = gg1 + reshape(be*X.node(k1,:),[L L fl]);
		gg2 = gg2 + reshape(be*X.node(k2,:),[L L fl]);
	  end;
	  
	  if C.symmetric
		gg1 = gg1 + permute(gg2,[2 1 3]);
		gg2(:) = 0;
	  end;
	  
	  G.dhEdge = G.dhEdge + I * cat(4,gg1,gg2);
	  
	  
	  % Vertical Features
	  szd = size(G.dhEdge);
	  
	  fl = size(X.node,2);
	  
	  gg1 = zeros(szd(1:3));
	  gg2 = zeros(szd(1:3));
	  
	  [R,M] = ind2sub([sz(1)-1 sz(2)],1:size(X.vEdge,1));
	
	  for e=1:length(R)
		
		k1 = K(R(e),M(e));
		k2 = K(R(e)+1,M(e));
		
        be = B.vEdge(e,:)'; % L^2x1

		
		gg1 = gg1 + reshape(be*X.node(k1,:),[L L fl]);
		gg2 = gg2 + reshape(be*X.node(k2,:),[L L fl]);
		
	  end;
	  
	  if C.symmetric
		gg1 = gg1 + permute(gg2,[2 1 3]);
		gg2(:) = 0;
	  end;
	  
	  G.dvEdge = G.dvEdge + I * cat(4,gg1,gg2);
	  
	  if C.isotropic
		G.dhEdge = G.dhEdge + G.dvEdge;
		G.dvEdge(:) = 0;
	  end;
	  
	end; % default edge features

  end; % edge params
  

	  
  
function [A,B] = order(M,N)
  A = min(M,N);
  B = max(M,N);
 
