function t = tauplus_sos(A)

% t = TAUPLUS_SOS(A)
%
% Computes the quantity \tau_+^{sos}(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 nonnegative matrix
%
% OUTPUT
%    t: a real number equal to \tau_+^{sos}(A)
%
% EXAMPLE
%    >> A = [1 1 0 0; 1 0 1 0; 0 1 0 1; 0 0 1 1];
%    >> tauplus_sos(A)   % Outputs 4
%
% NOTE
%    Requires Yalmip (http://users.isy.liu.se/johanl/yalmip/)
%
% AUTHOR
%    Hamza Fawzi (hfawzi@mit.edu)

[m,n] = size(A);

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

% Extract strictly positive elements of A
ix = find(A > 0);
a = A(ix);
N = numel(ix);

IX = zeros(size(A));
IX(ix) = 1:N;

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

% Setup the linear equalities X_{ij,kl} = X_{il,kj} corresponding
% to the rank-one constraint (vanishing of 2x2 minors)
for ii=1:m,
    for jj=1:n,
        for kk=(ii+1):m,
            for ll=(jj+1):n,
                ind_ij = IX(ii,jj);
                ind_kl = IX(kk,ll);
                ind_il = IX(ii,ll);
                ind_kj = IX(kk,jj);
                if A(ii,jj)*A(kk,ll) > 0 && A(ii,ll)*A(kk,jj) == 0
                    p_cons = [p_cons; X(ind_ij,ind_kl) == 0];
                elseif A(ii,jj)*A(kk,ll) == 0 && A(ii,ll)*A(kk,jj) > 0
                    p_cons = [p_cons; X(ind_il,ind_kj) == 0];
                elseif A(ii,jj)*A(kk,ll) > 0 && A(ii,ll)*A(kk,jj) > 0
                    p_cons = [p_cons; X(ind_ij,ind_kl) == X(ind_il,ind_kj)];
                end
            end
        end
    end
end

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

t = double(t);

end