function sinr_score = compute_sinr(A,X,B,z,e,Sigma)
    
    [~,k]          = size(X);
    s              = z'; % (k,n) dimensional
    eta            = e'; % (k,n) dimensional
    sinr_score     = 0;
    % sinr_opt_demix = A'*pinv(cov(X));
    sinr_opt_demix = A'*pinv(B*B' + Sigma);
    match_matrix   = zeros(k,k);

    for i=1:k
        for j=1:k
            u                 = sinr_opt_demix(j,:);
            num               = var(u*B(:,i)*s(i,:));
            denom             = var(u*B*s - u*B(:,i)*s(i,:)) + var(u*eta);
            score             = num/denom;
            match_matrix(i,j) = score;
        end
    end

    [~,cost] = hungarian_assignment(-match_matrix);
    sinr_score = -cost/k;
    % keyboard;
    %match_matrix
end