function Rr = restrict(ny,nk,ng,options)
% restrict  [Not a public function] Convert parameter restrictions to hyperparameter matrix form.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

%********************************************************************

if isempty(options.constraints) ...
      && isempty(options.a) ...
      && isempty(options.c) ...
      && isempty(options.g)
   Rr = [];
end

if isnumeric(options.constraints)
   Rr = options.constraints;
   return
end

nlag = options.order;
if options.diff
   nlag = nlag - 1;
end

nbeta = ny*(nk+ny*nlag+ng);
Q = zeros([0,nbeta]);
q = zeros(0);

isplain = ~isempty(options.a) ...
    || ~isempty(options.c) ...
    || ~isempty(options.g);

% General constraints.
rstring = '';
restrict = lower(strtrim(options.constraints));
if ~isempty(restrict)
   restrict = strfun.converteols(restrict);
   restrict = strrep(restrict,char(10),' ');
   restrict = lower(restrict);
   % Replace semicolons outside brackets with &s.
   restrict = strfun.strrepoutside(restrict,';','&','[]','()');
   % Read individual &-separated restrictions.
   eachrestrict = regexp(restrict,'(.*?)(?:&|$)','tokens');
   % Convert restrictions to implicit forms.
   eachrestrict = regexprep([eachrestrict{:}],...
      '=(.*)','-\($1\)');
   % Vectorise and vertically concatenate restrictions.
   for i = 1 : numel(eachrestrict)
      rstring = [rstring,'vec_(',eachrestrict{i},');'];
   end
   if ~isempty(rstring)
      rstring = ['[',rstring,']'];
   end
end

% A, C, G restrictions.
if ~isempty(rstring)
% General constraints exist. Set up (Q,q) first for general and plain
% constraints, then convert them to (R,r).
   rfunc = eval(['@(c,a,g) ',rstring,';']);
   [Q1,q1] = generalrestrict_(rfunc,ny,nk,ng,nlag);
   Q = [Q;Q1];
   q = [q;q1];
    % Plain constraints.
    if isplain
       [Q2,q2] = plainrestrict1_(options,ny,nk,ng,nlag);
       Q = [Q;Q2];
       q = [q;q2];
    end
    % Convert Q*beta + q = 0 to beta = R*gamma + r,
    % where gamma is a vector of free hyperparameters.
    if ~isempty(Q)
       R = null(Q);
       r = -pinv(Q)*q;
       Rr = sparse([R,r]);
    end
elseif isplain
   [R,r] = plainrestrict2_(options,ny,nk,ng,nlag);
   Rr = sparse([R,r]);
end

end


%********************************************************************
% Subfunction generalrestrict_().
function [Q,q] = generalrestrict_(rfunc,ny,nk,ng,nlag)
   % Q*beta = q
   aux = reshape(transpose(1:ny*(nk+ny*nlag+ng)),[ny,nk+ny*nlag+ng]);
   cindex = aux(:,1:nk);
   aux(:,1:nk) = [];
   aindex = reshape(aux(:,1:ny*nlag),[ny,ny,nlag]);
   aux(:,1:ny*nlag) = [];
   gindex = aux;
   c = zeros(size(cindex)); % Constant.
   a = zeros(size(aindex)); % Transition matrix.
   g = zeros(size(gindex)); % Cointegrating vector.
   % Q*beta + q = 0.
   try
      q = rfunc(c,a,g);
   catch Error
      VAR.error(1,{Error.message});
   end
   nrestrict = size(q,1);
   Q = zeros([nrestrict,ny*(nk+ny*nlag+ng)]);
   for i = 1 : numel(c)
      c(i) = 1;
      Q(:,cindex(i)) = rfunc(c,a,g) - q;
      c(i) = 0;
   end
   for i = 1 : numel(a)
      a(i) = 1;
      Q(:,aindex(i)) = rfunc(c,a,g) - q;
      a(i) = 0;
   end
   for i = 1 : numel(g)
       g(i) = 1;
       Q(:,gindex(i)) = rfunc(c,a,g) - q;
       g(i) = 0;
   end
end
% End of subfunction generalrestrict_().

%********************************************************************
% Subfunction plainrestrict1_().
function [Q,q] = plainrestrict1_(options,ny,nk,ng,nlag)
   [A,C,G] = assignplainrestrict_(options,ny,nk,ng,nlag);
   nbeta = ny*(nk+ny*nlag+ng);
   % Construct parameter restrictions first,
   % Q*beta + q = 0,
   % splice them with the general restrictions
   % and only then convert these to hyperparameter form.
   Q = eye(nbeta);
   q = -[C,A(:,:),G];
   q = q(:);
   index = ~isnan(q);
   Q = Q(index,:);
   q = q(index);
end
% End of subfunction plainrestrict1_().

%********************************************************************
% Subfunction plainrestrict2_().
function [R,r] = plainrestrict2_(options,ny,nk,ng,nlag)
   [A,C,G] = assignplainrestrict_(options,ny,nk,ng,nlag);
   nbeta = ny*(nk+ny*nlag+ng);
   % Construct directly hyperparameter form:
   % beta = R*gamma + r.
   R = eye(nbeta);
   r = [C,A(:,:),G];
   r = r(:);
   index = ~isnan(r);
   R(:,index) = [];
   r(~index) = 0;
end
% End of subfunction plainrestrict2_().

%********************************************************************
% Subfunction assignplainrestrict_().
function [A,C,G] = assignplainrestrict_(options,ny,nk,ng,nlag)
   A = nan([ny,ny,nlag]);
   C = nan([ny,nk]);
   G = nan([ny,ng]);
   if ~isempty(options.a)
      try
         A(:,:,:) = options.a;
      catch
         VAR.error(35,{'A',sprintf('%g-by-%g-by-%g',ny,ny,nlag)});
      end
   end
   if ~isempty(options.c)
      try
         C(:,:) = options.c;
      catch
         VAR.error(35,{'C',sprintf('%g-by-%g',ny,nk)});
      end
   end
   if ~isempty(options.g)
      try
         G(:,:) = options.g;
      catch
         VAR.error(35,{'G',sprintf('%g-by-%g',ny,ng)});
      end
   end
end
% End of subfunction assignplainrestrict_().

%********************************************************************
% Subfunction vec_().
function x = vec_(x)
   x = x(:);
end
% End of subfunction vec_().