function rho = rho_calc (rho_output, leave_one_out)
%RHO_CALC  Convert output from nnrho.cpp to rho values
%
%   rho = RHO_CALC(rho_output);

if nargin < 2, leave_one_out = 1; end

rho_output = sum(rho_output, 1);
nx = (size(rho_output,2) - 1) / 3;
n = sum(rho_output(nx+nx+1:end));
disp([num2str(nx) ' samples, ' num2str(n) ' trials']);

mean_ei   = rho_output(1:nx) / n;
% we omit 1/nx below since a constant doesn't affect corrcoef
mean_eipi = rho_output(nx+1:nx+nx) / n;
mean_pi   = sum(mean_ei);
mean_pi2  = sum(mean_eipi);
mean_ei2  = mean_ei;

if leave_one_out,
    mean_pi2 = mean_pi2 - 2*mean_eipi + mean_ei2;
    mean_eipi = mean_eipi - mean_ei2;
    mean_pi = mean_pi - mean_ei;
end

sigma_ei = sqrt(mean_ei2 - mean_ei.^2);
sigma_pi = sqrt(mean_pi2 - mean_pi.^2);

rho = (mean_eipi - mean_ei.*mean_pi) ./ (sigma_ei.*sigma_pi);