function [out] = shap_loss_comp(x_target,X_bg,f,mod_params,loss_params,id_const,bin_edges,method,comp_p0,TVP,verbose)
    

    [n_bg] = size(X_bg,1); % number of background samples

    B=3;
    bins_center = 0.5 * (bin_edges(1:end-1) + bin_edges(2:end));

    d      = size(X_bg,2);
    vars   = 1:d;
    if isempty(id_const)
        d      = size(X_bg,2);    % number of variables
    else
        % constants have to be ordered first
        d      = sum(~ismember(vars,id_const));  % new variable count w/o constants
        ids    = vars(~ismember(vars,id_const)); % remove constants given 0 contributions
    end
    
    % construct permutations
    P = perms(1:d)+length(id_const);

    % construct baseline
    if comp_p0 == 1 || isempty(comp_p0)
        p_bg = zeros(n_bg, 3);
        if  TVP == 0 || isempty(TVP)
            for j = 1:n_bg
                [p_bg(j, :)] = f(X_bg(j,:),mod_params,loss_params);
            end
        else
            for j = 1:n_bg
                [p_bg(j, :)] = f(X_bg(j,:),mod_params(end,:),loss_params);
            end
            mod_params = mod_params(end,:);
        end
        phi0 = mean(p_bg, 1);
        out.phi0 = phi0;
    end

    % construct forecast
    p_full     = f(x_target,mod_params,loss_params);
    out.p_full = p_full;

    phi = zeros(d+length(id_const), 3);
    if method == "exact"
        for p = 1:size(P,1)
            order = P(p, :);
            for j = 1:n_bg
                xj = X_bg(j, :);
                xS = xj;
                for i = 1:d
                    feat = order(i);
                    xS_feat = xS;
                    xS_feat(feat) = x_target(feat);
                    p1 = f(xS,mod_params,loss_params);
                    p2 = f(xS_feat,mod_params,loss_params);
                    phi(feat, :) = phi(feat, :) + (p2 - p1) / (size(P,1) * n_bg);
                    xS = xS_feat;
                end
            end
        end
        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
    elseif method == "random"
        M = n_bg*200;1000;
        for m = 1:M
            if isempty(id_const)
                order = randsample(vars,d);%randperm(d);
            else
                order = randsample(vars(vars~=id_const),d);%randperm(d);
            end
            j = randi(n_bg);
            xj = X_bg(j, :);
            xS = xj;
            for i = 1:d
                feat = order(i);
                xS_feat = xS;
                xS_feat(feat) = x_target(feat);
                p1 = f(xS,mod_params,loss_params);
                p2 = f(xS_feat,mod_params,loss_params);
                phi(feat, :) = phi(feat, :) + (p2 - p1) / M;
                xS = xS_feat;
            end
        end
        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
    elseif method == "kernel"
         M = 1000;
        Z = [];
        Y = [];
        W = [];

        if isempty(id_const)
            id_const = 0;
        end

        for m = 1:M
            S = [zeros(1,id_const) randi([0 1], 1, d)];%randi([0 1], 1, d);
            if all(S==0) || all(S==1)
                continue;
            end

            j = randi(n_bg);
            xj = X_bg(j, :);
            x_S = xj;
            x_S(logical(S)) = x_target(logical(S));

            for i = 1+id_const:d+id_const
                if ~S(i)
                    x_Si = x_S;
                    x_Si(i) = x_target(i);
                    delta  = f(x_Si,mod_params,loss_params) - f(x_S,mod_params,loss_params);
                     % Store subset indicator, delta, and weight
                    z_row = zeros(1, d);
                    z_row(i) = 1;
                    Z = [Z; z_row];
                    Y = [Y; delta];

                    s = sum(S);
                    w = factorial(s) * factorial(d - s - 1) / factorial(d);
                    W = [W; w];
                end
            end
        end

        for b = 1:B
            y_b = Y(:, b);
            phi_b = (Z' * diag(W) * Z) \ (Z' * diag(W) * y_b);
            phi(:, b) = phi_b;
        end

        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
    elseif method == "knn"
        K = round(0.1*size(X_bg,1));
        M = 1000;
        distances = vecnorm(X_bg - x_target, 2, 2);
        [~, idxs] = mink(distances, K);
        X_knn = X_bg(idxs, :);
        for m = 1:M
            if isempty(id_const)
                order = randsample(vars,d);%randperm(d);
            else
                order = randsample(vars(vars~=id_const),d);%randperm(d);
            end
            xj = X_knn(randi(K), :);
            xS = xj;
            for i = 1:d
                feat = order(i);
                xS_feat = xS;
                xS_feat(feat) = x_target(feat);
                p1 = f(xS,mod_params,loss_params);
                p2 = f(xS_feat,mod_params,loss_params);
                phi(feat, :) = phi(feat, :) + (p2 - p1) / M;
                xS = xS_feat;
            end
        end
        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
        
        for i = 1:size(X_knn,1)
            aux(i,:) = f(X_knn(i,:),mod_params,loss_params);
        end

        out.phi0cond = mean(aux, 1);  % Lo
    end

    % Plot forecast and background density
    % if isempty(verbose) || verbose == 0
    %     figure;
    %     subplot(2,1,1);
    %     plot(bins_center, p_full, 'b-', 'LineWidth', 2); hold on;
    %     plot(bins_center, mean(p_bg), 'k--', 'LineWidth', 2);
    %     title('Forecast vs Background Densities');
    %     legend('Forecast at target x', 'Background avg density');
    %     grid on;
    % 
    %     subplot(2,1,2);
    %     bar_data = phi';
    %     bar_data = bar_data(:,~sum(bar_data)==0);
    %     bar(bins_center, bar_data, 'stacked');
    %     hold on;
    %     plot(bins_center, p_full - phi0, 'ko--', 'LineWidth', 1.5);
    %     title('SHAP values for bin probabilities');
    %     xlabel('Bin Center'); ylabel('Δ Probability');
    %     legend(arrayfun(@(i) sprintf('Feature %d', i), 1:d, 'UniformOutput', false),'Location','NorthWest');
    %     grid on;
    % 
    % end

end