import numpy as np


def permutation_importance(
    model_predict, X, y, metric_fn, n_repeats=10, seed=0
):
    rng = np.random.default_rng(seed)
    baseline_score = metric_fn(y, model_predict(X))
    importances = np.zeros(X.shape[1])

    for feature_idx in range(X.shape[1]):
        drops = []
        for _ in range(n_repeats):
            X_permuted = X.copy()
            X_permuted[:, feature_idx] = rng.permutation(
                X_permuted[:, feature_idx]
            )
            permuted_score = metric_fn(y, model_predict(X_permuted))
            drops.append(baseline_score - permuted_score)
        importances[feature_idx] = np.mean(drops)

    return importances


def kernel_shap_approx(model_predict, x, background, n_samples=500, seed=0):
    """A simplified KernelSHAP: approximate Shapley values via
    weighted linear regression over random feature-subset masks."""
    rng = np.random.default_rng(seed)
    n_features = len(x)
    baseline = background.mean(axis=0)

    masks, outputs = [], []
    for _ in range(n_samples):
        mask = rng.integers(0, 2, size=n_features)
        perturbed = np.where(mask == 1, x, baseline)
        masks.append(mask)
        outputs.append(model_predict(perturbed.reshape(1, -1))[0])

    masks = np.array(masks, dtype=float)
    outputs = np.array(outputs)

    # simplified SHAP kernel weighting (Shapley kernel)
    def shap_kernel_weight(mask):
        z = mask.sum()
        if z == 0 or z == n_features:
            return 1e6  # empty/full coalitions anchor the regression
        from math import comb
        return (n_features - 1) / (
            comb(n_features, int(z)) * z * (n_features - z)
        )

    weights = np.array([shap_kernel_weight(m) for m in masks])

    # weighted least squares: masks @ phi ≈ outputs - f(baseline)
    f_baseline = model_predict(baseline.reshape(1, -1))[0]
    target = outputs - f_baseline
    W = np.diag(weights)
    A = masks.T @ W @ masks + 1e-6 * np.eye(n_features)
    b = masks.T @ W @ target
    phi = np.linalg.solve(A, b)
    return phi


def lime_explain(model_predict, x, n_samples=500, kernel_width=0.75, seed=0):
    rng = np.random.default_rng(seed)
    n_features = len(x)
    perturbations = rng.normal(0, 1, size=(n_samples, n_features))
    perturbed_x = x + perturbations * 0.3
    predictions = model_predict(perturbed_x)

    distances = np.linalg.norm(perturbations, axis=1)
    weights = np.exp(-(distances ** 2) / (kernel_width ** 2))

    # weighted linear regression: predictions ~ perturbations
    W = np.diag(weights)
    A = perturbations.T @ W @ perturbations + 1e-6 * np.eye(n_features)
    b = perturbations.T @ W @ predictions
    coefficients = np.linalg.solve(A, b)
    return coefficients
