function t = tauplus_sos_tensor(A)

% t = TAUPLUS_SOS_TENSOR(A)
%
% Computes the quantity \tau_+^{sos}(A) for a 3-dimensional tensor A as
% defined in the paper
%
%   "Self-scaled bounds for atomic cone ranks:
%   applications to nonnegative rank and cp-rank"
%   Hamza Fawzi, Pablo A. Parrilo
%   http://arxiv.org/abs/1404.3240
%
% This is a lower bound on the nonnegative rank of A.
%
% INPUT
%    A: a 3-dimensional nonnegative array
%
% OUTPUT
%    t: a real number equal to \tau_+^{sos}(A)
%
% EXAMPLE
%    >> A = zeros(2,2,2); A(:,:,1) = [1 1; 1 0]; A(:,:,2) = [0 1; 1 1];
%    >> tauplus_sos_tensor(A)   % Outputs 3
%
% NOTE
%    Requires Yalmip (http://users.isy.liu.se/johanl/yalmip/)
%
% AUTHOR
%    Hamza Fawzi (hfawzi@mit.edu)

if any(A(:) < 0)
	error('A has to be nonnegative');
end

if numel(size(A)) ~= 3,
    error('A has to be a 3-dimensional array');
end

a = A(:);
[n1,n2,n3] = size(A);
N = n1*n1*n3;

X = sdpvar(N,N,'symmetric');
t = sdpvar(1);
p_cons = [[t a'; a X] >= 0; diag(X) <= a.^2];

% Setup the linear equalities corresponding to the rank-one conditions for
% 3-dimensional tensors
s2i = @(x,y,z) sub2ind(size(A),x,y,z);
if n1 == 2 && n2 == 2 && n3 == 2,
    % Use hard-coded equations for the 2x2x2 case (9 equations in total)
    p_cons = [
        p_cons;
        % 2x2 minor in slice (:,:,1)
        X(s2i(1,1,1),s2i(2,2,1)) == X(s2i(1,2,1),s2i(2,1,1));
        % 2x2 minor in slice (:,:,2)
        X(s2i(1,1,2),s2i(2,2,2)) == X(s2i(1,2,2),s2i(2,1,2));
        % 2x2 minor in slice (:,1,:)
        X(s2i(1,1,1),s2i(2,1,2)) == X(s2i(1,1,2),s2i(2,1,1));
        % 2x2 minor in slice (:,2,:)
        X(s2i(1,2,1),s2i(2,2,2)) == X(s2i(1,2,2),s2i(2,2,1));
        % 2x2 minor in slice (1,:,:)
        X(s2i(1,1,1),s2i(1,2,2)) == X(s2i(1,1,2),s2i(1,2,1));
        % 2x2 minor in slice (2,:,:)
        X(s2i(2,1,1),s2i(2,2,2)) == X(s2i(2,1,2),s2i(2,2,1));
        % 2x2 minors between (1,1,1) and (2,2,2)
        X(s2i(1,1,1),s2i(2,2,2)) == X(s2i(1,1,2),s2i(2,2,1));
        X(s2i(1,1,1),s2i(2,2,2)) == X(s2i(1,2,1),s2i(2,1,2));
        X(s2i(1,1,1),s2i(2,2,2)) == X(s2i(2,1,1),s2i(1,2,2));
    ];
else
    % Generate equations in the general case
    % Note: this code is not optimized and there may be redundant equations
    lst = [0 0 0 0];
    % Loop over multi-indices alpha and beta
    for ind_alpha=1:N,
        for ind_beta=(ind_alpha+1):N,
            alpha = cell(1,3);
            beta = cell(1,3);
            [alpha{:}] = ind2sub(size(A),ind_alpha);
            [beta{:}] = ind2sub(size(A),ind_beta);
            alpha = [alpha{:}];
            beta = [beta{:}];
            for k=1:3,
                if alpha(k) ~= beta(k)
                    % Swap index k of alpha and beta
                    alpha1 = alpha; alpha1(k) = beta(k);
                    beta1 = beta; beta1(k) = alpha(k);
                    ind_alpha1 = s2i(alpha1(1),alpha1(2),alpha1(3));
                    ind_beta1 = s2i(beta1(1),beta1(2),beta1(3));
                    sig = [ind_alpha ind_beta ...
                            min(ind_alpha1,ind_beta1) ...
                            max(ind_alpha1,ind_beta1)];
                    if ~ismember(sig,lst,'rows')
                        lst = union(lst,sig,'rows');
                        % Add constraint
                        p_cons = [p_cons; X(ind_alpha,ind_beta) == ...
                            X(ind_alpha1,ind_beta1)];
                    end
                end
            end
        end
    end
    fprintf('# of equations added = %d\n', size(lst,1)-1);
end

p_obj = t;
p_sol = solvesdp(p_cons,p_obj);

t = double(t);

end