function G = unrollsubgm(W,N,H,V,I,D)
% UNROLLSUBGM Unroll a grid CRF for data creating a graphical sub-model
%
% G = UNROLLSUBGM(C,N,H,V,I) where C is the gcrf object N,H,V are the
% node, horizontal and vertical feature vectors, and I is the list of
% nodes that are included in the model.

            
% Jerod Weinman
% jerod@acm.org
  
%	$Id: unrollsubgm.m 1 2009-05-07 13:41:00Z 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/>.

	
  if nargin<6
	D = 1;
  end;
  
  sz = gridsize(N,H);

  L = size(W.node,2);

  if size(I,1)~=length(I)
	I = I';
  end;
  
  % Node inclusion indicator
  A = sparse(size(N,1),1);
  A(I) = 1:length(I);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Allocate adjacency matrix

  K = nodenum(sz(1), sz(2));
  
  E = logical(sparse(length(I),length(I)));
  

  S = cell(size(E));
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate node energies

  % Note: notes are ordered in the grid with their standard matlab indices
  
  % Perform standard multiplication (weight dot prod)
  P = N(I,:)*W.node(1:end-1,:) + ...
	  ones(length(I),1)*W.node(end,:);
  
  % Transpose for inferencer
  P = P';
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate horizontal energies

  % H:  NxD
  % We: DxL^2
  % Uh: NxL^2
  We = reshape(W.hEdge,[L*L size(W.hEdge,3)])';


  % Only include edges that are connected to two included nodes
  
  E0 = mkrvector(K(:,1:end-1)); % left node list
  E1 = mkrvector(K(:,2:end)); % right node list
  
  Ie = find(A(E0) & A(E1));
      
  % Weight dot product
  Uh = H(Ie,:)*We(1:end-1,:) + ones(length(Ie),1)*We(end,:);

  if D
	
	% Note: We could reduce the computation here by checking to see
    % if the CRF is symmetric and/or isotropic
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Calculate default horizontal energies
	
	% Node Data: NxD
	% Wex: DxL^2
	% Un: NxL^2
	We1 = reshape(W.dhEdge(:,:,:,1), ...
				  [L*L size(W.dhEdge,3)])';
	
	We2 = reshape(W.dhEdge(:,:,:,2), ...
				  [L*L size(W.dhEdge,3)])';

	% Weight dot product
	Un1 = N(I,:)*We1;
	Un2 = N(I,:)*We2;
  end;

  
  for e=1:length(Ie)
	k1 = A(E0(Ie(e)));
	k2 = A(E1(Ie(e)));

	E(k1,k2) = 1;
	
	psi = reshape(Uh(e,:)',[L L]);
	
	if D
	  psi = psi + reshape((Un1(k1,:)+Un2(k2,:))',[L L]);
	end;
	
	S{k1,k2} = psi;
	
  end;
  


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate vertical energies

  % V: NxD
  % We: DxL^2
  % Uv: NxL^2

  We = reshape(W.vEdge,[L*L size(W.vEdge,3)])';
	
  % Only include edges that are connected to two included nodes
  
  E0 = mkrvector(K(1:end-1,:));
  E1 = mkrvector(K(2:end,:));
  
  Ie = find(A(E0) & A(E1));
  
  % Weight dot product
  Uv = V(Ie,:)*We(1:end-1,:) + ones(length(Ie),1)*We(end,:);

  if D
		
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Calculate default vertical energies

	% Node Data: NxD
	% Wex: DxL^2
	% Un: NxL^2
	We1 = reshape(W.dvEdge(:,:,:,1), ...
				  [L*L size(W.dvEdge,3)])';
	
	We2 = reshape(W.dvEdge(:,:,:,2), ...
				  [L*L size(W.dvEdge,3)])';

	% Weight dot product
	Un1 = N(I,:)*We1;
	Un2 = N(I,:)*We2;
  end;

  for e=1:length(Ie)

	k1 = A(E0(Ie(e)));
	k2 = A(E1(Ie(e)));

	E(k1,k2) = 1;
		
	psi = reshape(Uv(e,:)',[L L]);
	
	if D
	  psi = psi + reshape((Un1(k1,:)+Un2(k2,:))',[L L]);
	end;

	S{k1,k2} = psi;
	
  end;
  
  
  G = struct('numNodes', {length(I)}, ...
			 'adjmat', {E}, ...
			 'numLabels', {L}, ...
			 'nodeEnergies',{P}, ...
			 'edgeEnergies',{S} );
