function t = taucp_sos(A)

% t = TAUCP_SOS(A)
%
% Computes the quantity \tau_{cp}^{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 cp-rank of A.
%
% INPUT
%    A: a completely positive matrix
%
% OUTPUT
%    t: a real number equal to \tau_{cp}^{sos}(A)
%
% EXAMPLE
%    >> A = [3 0 1 1 1; 0 3 1 1 1; 1 1 2 0 0; 1 1 0 2 0; 1 1 0 0 2];
%    >> taucp_sos(A)   % Outputs 6
%
% NOTE
%    Requires Yalmip (http://users.isy.liu.se/johanl/yalmip/)
%
% AUTHOR
%    Hamza Fawzi (hfawzi@mit.edu)


if size(A,1) ~= size(A,2)
	error('A has to be a symmetric matrix');
end

n = size(A,1);

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

a = A(:);
N = length(a);

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

% Setup the linear equalities X_{ij,kl} = X_{il,kj} corresponding
% to the rank-one constraint (vanishing of 2x2 minors)
for ii=1:n,
    for jj=1:n,
        for kk=(ii+1):n,
            for ll=(jj+1):n,
                ind_ij = sub2ind([n n],ii,jj);
                ind_kl = sub2ind([n n],kk,ll);
                ind_il = sub2ind([n n],ii,ll);
                ind_kj = sub2ind([n n],kk,jj);
                p_cons = [p_cons; X(ind_ij,ind_kl) == X(ind_il,ind_kj)];
            end
        end
    end
end
p_obj = t;
p_sol = solvesdp(p_cons,p_obj);

t = double(t);

end