function [ kappa, A, T ] = cpdcond( A, ~ )
%CPDCOND Computes the condition number. 
%  The routine estimates the norm-balanced condition number of the factor
%   matrices A. If the second input argument is specified, the absolute
%   condition number is computed. 
%  The first output argument is the norm-balanced condition number.
%  The second output argument are the norm-balanced factor matrices with
%   respect to which the condition number is computed.
%  The third output argument contains Terracini's matrix of the
%   norm-balanced representatives.

% Constants
d = length(A);
n = arrayfun(@(k) size(A{k},1), 1:d);
r = size(A{1},2);
Pi = prod(n);
Sigma = sum(n-1);

%% Apply norm-balancing
scal = ones(1,r);
for k = 1 : d
    lscal = arrayfun(@(l) norm(A{k}(:,l)), 1:r);
    scal = scal .* lscal.^(1/d);
    A{k} = A{k} * diag(1./lscal);
end
for k = 1 : d
    A{k} = A{k} * diag(scal);
end
normP = sqrt(d)*norm(scal);

% Condition number is unbounded by definition.
if r*(1+Sigma) > Pi
   warning('The number of terms is not subgeneric.')
   kappa = inf;
   T = [];
   return;
end

%% Construct Terracini's matrix
T = zeros(Pi,r*(d+Sigma));
for l = 1 : r
   lT = [];
   for k = 1 : d
       B = arrayfun(...
           @(kk) repmat(A{kk}(:,l),1,n(k)), 1:d, 'UniformOutput', false ...
       );
       B{k} = eye(n(k));
       lT = [lT krp(B)];
   end
   T(:,(d+Sigma)*(l-1)+1:(d+Sigma)*l) = lT;
end

%% Compute relative condition number
N = r*(1+Sigma);
S = svd(T,0);
kappa = 1 / S(N);
if nargin == 1
   tensNorm = norm(krp(A)*ones(r,1));
   kappa = (kappa * tensNorm) / normP;
end
if ( S(1)*100*eps > S(N) )
   warning('Computed condition number may not be accurate!')
end

    % Compute Khatri-Rao product
    function [X] = krp(B)
       Ns = prod(arrayfun(@(k)size(B{k},1),1:length(B)));
       X = zeros(Ns,size(B{1},2));
       for ll = 1 : size(X,2)
           x = B{1}(:,ll);
           for kk = 2 : length(B)
               x = kron(x,B{kk}(:,ll));
           end
           X(:,ll) = x;
       end
    end
end

