function [F,error,R2] = RegrMLRC(X,Y,G,g,theta)

%   [F,error,R2] = RegrMLRC(X,Y,G,g,theta)
%   [F,error,R2] = RegrMLRC(X,Y,G,g)
%
% Regression from X to Y with constraint G*F = g.
% Each column in F has the same constraints, different g.
%
% Input parameters:
%  - X: Input data block (k x n)
%  - Yi: Output data block (k x m)
%  - G,g: Linear constraints in the form GF = g
%  - theta: Latent basis, orthogonal or not (optional)
% Return parameters:
%  - F: Mapping matrix, Yhat = X*F
%  - error: Prediction errors
%  - R2: Data fitting criterion R^2
%
% Heikki Hyotyniemi Dec.21, 2002


F = NaN; error = NaN; R2 = NaN; 
[kx,n] = size(X);
[ky,m] = size(Y);
if kx ~= ky, disp('Incompatible X and Y'); break; 
else k = kx; end
if size(G,1) ~= size(g,1), disp('Incompatible G and g'); break; end
if size(G,2) ~= n, disp('Incompatible G and X'); break; end
if size(g,2) ~= m, disp('Incompatible g and Y'); break; end

if nargin<5
   invXX = inv(X'*X);
   F = invXX*(X'-G'*inv(G*invXX*G')*G*invXX*X')*Y ...
     + invXX*G'*inv(G*invXX*G')*g;
else
   F1 = theta*inv(theta'*theta);
   Z = X*F1;
   invZZ = inv(Z'*Z);
   F2 = invZZ*(Z'-F1'*G'*inv(G*F1*invZZ*F1'*G')*G*F1*invZZ*Z')*Y ...
        + invZZ*F1'*G'*inv(G*F1*invZZ*F1'*G')*g;
   F = F1*F2;
end

Yhat = X*F;
error = Yhat - Y;
R2 = 1 - sum((error).^2)./sum(Y.^2);