function [fid,X,D,N] = fid_sio(rho)

d = size(rho,2);
if d==1
    rho=rho*rho';
    d=size(rho,1);
end
rho=(rho+rho')/2;

cvx_begin sdp quiet

variable X(d,d) hermitian
dual variable D
dual variable N

X>= -eye(d)
X<= eye(d)
D: diag(X) == 0
N: X(:) == nonnegative(d^2)

maximize (trace(abs(rho)*X))

cvx_end

fid = (cvx_optval+1)/2;
D = sparse(diag(D));
N = reshape(N, [d,d]);
N = (N + N')/2;

end