function [Xr, Sr, rho, eta] = pttls(V, d, missing, r, Sglobal)
%PTTLS Truncated TLS regularization with permuted columns.
%
%    Given matrices A and B, the total least squares (TLS) problem
%    consists of finding a matrix Xr that satisfies
%
%                (A+dA)*Xr = B+dB.
%
%    The solution must be such that the perturbation matrices dA
%    and dB have minimum Frobenius norm rho=norm( [dA dB], 'fro')
%    and each column of B+dB is in the range of A+dA [1].
%  
%    [Xr, Sr, rho, eta] = PTTLS(V, d, colA, colB, r) computes the
%    minimum-norm solution Xr of the TLS problem, truncated at rank r
%    [2]. The solution Xr of this truncated TLS problem is a
%    regularized error-in-variables estimate of regression
%    coefficients in the regression model A*X = B + noise(S). The
%    model may have multiple right-hand sides, represented as columns
%    of B.
%
%    As input, PTTLS requires the right singular matrix V of the
%    augmented data matrix C = U*diag(s)*V' and the vector d=s.^2 with
%    the squared singular values. Only right singular vectors V(:,j)
%    belonging to nonzero singular values s(j) are required.  Usually,
%    the first n columns of the augmented data matrix C correspond to
%    the n columns of the matrix A, and the k last columns to the k
%    right-hand sides B, so that the augmented data matrix is of the
%    form C=[A B]. PTTLS allows a more flexible composition of the
%    data matrix C: the columns with indices colA correspond to
%    columns of A; the columns with indices colB correspond to columns
%    of B.
%
%    The right singular vectors V and the squared singular values d
%    may be obtained from an eigendecomposition of the matrix C'*C = v
%    * d * v', which, for centered data, is proportional to the sample
%    covariance matrix.
%
%    PTTLS returns the rank-r truncated TLS solution Xr and the matrix
%    Sr = dB'*dB, which is proportional to the estimated covariance
%    matrix of the residual dB. Also returned are the Frobenius norm 
%
%                    rho = norm([dA dB], 'fro') 
%
%    of the residuals and the Frobenius norm 
%
%                    eta = norm(Xr,'fro') 
%
%    of the solution matrix Xr.
%  
%    If the truncation parameter r(1:nr) is a vector of length nr,
%    then Xr is a 3-D matrix with
%    
%         Xr(:,:, 1:nr) = [ Xr(r(1)), Xr(r(2)), ..., Xr(r(nr)) ] .
%
%    The covariance matrix estimate Sr has an analogous structure, and
%    the residual norm rho(1:nr) and the solution norm eta(1:nr) are
%    vectors with one element for each r(1:nr).
%
%    If r is not specified or if r > n, r = n is used.

%     References: 
%     [1] Van Huffel, S. and J.Vandewalle, 1991:
%         The Total Least Squares Problem: Computational Aspects
%         and Analysis. Frontiers in Mathematics Series, vol. 9. SIAM.
%     [2] Fierro, R. D., G. H. Golub, P. C. Hansen and D. P. O'Leary, 1997:
%         Regularization by truncated total least squares, SIAM
%         J. Sci. Comput., 18, 1223-1241
%     [3] Golub, G. H, and C. F. van Loan, 1989: Matrix
%         Computations, 2d ed., Johns Hopkins University Press,
%         chapter 12.3 

error(nargchk(2,5,nargin))          % check number of input arguments 

[na,ma]   = size(V); 
nd        = length(d);
% if nd ~= prod(size(d))
%     error('Squared singular values d must be given as vector.')
% end
if ma < nd
    error(['All right singular vectors with nonzero singular value' ...
    ' are required.'])
end
d         = d(:);                   % make sure d is column vector
  
if nargin == 2
    n       = na-1;	              % default number of columns of A
    k       = 1;                      % default number of right-hand sides
    colA    = 1:n;                    % take first n columns of [A B] as A 
    colB    = na;                     % take last column of [A B] as B
else
    k       = nnz(missing);           % number of right-hand sides 
    n       = na-k;           % number of columns of A (number of variables)
end	     

if (nargin < 5) 
    r = n;                            % default truncation of TLS 
end

if (numel(r) ~= 1)
    error('Exactly one truncation parameter required.');
end


%   nr = length(r);                     % number of truncation parameters
%   if any(r < 1)
%     error('Impossible truncation parameter.')
%   end
  
% ir = find(r > n);
% if ~isempty(ir)
%     r(ir) = n;
%     warning('Truncation parameter lowered.')
% end

% initialize output variables
if (nargout > 2)
    rho = zeros(nr,1);
end

if (nnz(~missing) < r)
    r = max(min(r,n),1);
    Sr = V(missing,r+1:end)*(repmat(d(r+1:end),1,k).*V(missing,r+1:end)');
else    
    Sr = Sglobal(missing, missing);
end


V11 = V(~missing, 1:r);
V21 = V(missing, 1:r);
Xr = (V11 / (V11'*V11)) * V21';
    
% estimated covariance matrix of residuals dB0'*dB0 (up to a
% scaling factor)
%     V22         = V(colB, rc+1:nd); 
%     Sr(:,:,ir)  = V22 * (repmat(d(rc+1:nd), 1, k) .* V22');

%     Sr(:,:,ir)  =  V22 * diag(d(rc+1:end))* V22';
    
    
    % alternative formula when all left singular vectors are given:
    % V12         = V(colA, rc+1:end);
    % Xr(:,:,ir)  = -V12/V22; 

if (nargout > 2)       % residual norm = norm( [dA0 dB0], 'fro')
    rho   = sqrt(sum(d(r+1:nd)));     
end

if (nargout == 4)      % solution norm = norm( Xr, 'fro')
    eta   = norm(Xr, 'fro'); 
end