Merge pull request #32 from HashirA123/Probe

Probe Implementation

Author: zkhotanlou

experiments/experimental_setup.yaml

@@ -138,5 +138,7 @@ recourse_methods:
     hyperparams:
       loss_type: "BCE"
       binary_cat_features: True
+  probe:
+    hyperparams:
   roar:
     hyperparams:

experiments/results.csv

@@ -1415,6 +1415,34 @@ claproar,twomoon,linear,0.0,2.2577137270829443e-08,3.667492287085482e-16,1.87680
 claproar,twomoon,linear,0.0,1.6213392628472434e-08,1.431988503656087e-16,1.0531753025233572e-08,0.0,0.0,,,
 claproar,twomoon,linear,0.0,3.312241828035134e-08,7.202236947752843e-16,2.5826099814274528e-08,0.0,0.0,,,
 claproar,twomoon,linear,0.0,3.2146713291325575e-08,6.212848124877535e-16,2.330451231991049e-08,0.0,0.0,,,
+probe,adult,linear,51.0,2.5574682458217044,0.20724095760647737,0.10291039943695068,2.0,51.0,0.0,1.0,11.27795516
+probe,adult,linear,48.0,1.6209024338863478,0.0782517097074313,0.06244194507598877,2.0,51.0,,,
+probe,adult,linear,51.0,6.151970284269187,1.3522899002686382,0.268756240606308,2.0,48.0,,,
+probe,adult,linear,47.0,4.1338594518437715,0.6257325925713797,0.1863815188407898,2.0,44.0,,,
+probe,adult,linear,49.0,6.529014715911816,1.5907548243066623,0.2905552387237549,2.0,27.0,,,
+probe,compass,linear,7.0,1.0493243297463968,0.15758896097912403,0.15593880414962769,5.0,3.0,0.0,1.0,3.9818648800000007
+probe,compass,linear,7.0,1.1968591064214706,0.20529106263232566,0.178464874625206,5.0,3.0,,,
+probe,compass,linear,7.0,0.5305286708631014,0.040413185530297796,0.0809311717748642,5.0,6.0,,,
+probe,compass,linear,7.0,2.068292945623398,0.6289491666450724,0.3217114806175232,5.0,2.0,,,
+probe,compass,linear,6.0,0.11717507370600574,0.002380937481501453,0.024646831676363945,5.0,7.0,,,
+probe,credit,linear,20.0,0.7460355385040021,0.03599241548631901,0.05706929787993431,5.0,18.0,0.2666666666666667,1.0,4.403598066666665
+probe,credit,linear,18.0,0.2715397661723585,0.005983706032008009,0.028109369799494743,5.0,20.0,,,
+probe,credit,linear,18.0,0.2347176402264215,0.004717934816094677,0.025857295840978622,5.0,20.0,,,
+probe,german,linear,4.0,0.34145124321402176,0.02916159435279554,0.08767480622319615,2.0,2.0,0.19999999999999996,1.0,3.3755328400000026
+probe,german,linear,4.0,0.3343678011115066,0.027964667531387127,0.08589213144253283,2.0,2.0,,,
+probe,german,linear,4.0,0.267133860467632,0.017855097803688207,0.06914474771303289,2.0,2.0,,,
+probe,german,linear,4.0,0.3044240432907651,0.023182781515892512,0.07841214362312765,2.0,2.0,,,
+probe,german,linear,4.0,0.2473267066867472,0.015308234751386422,0.06424275391242085,2.0,3.0,,,
+probe,mortgage,linear,2.0,2.596893806724318,3.3737499230649304,1.3286230641796841,0.0,0.0,0.0,1.0,7.078433219999999
+probe,mortgage,linear,2.0,2.7662955305674433,3.8262065210593192,1.3854972163486525,0.0,0.0,,,
+probe,mortgage,linear,2.0,2.9121862424407814,4.241009224725797,1.473339416403974,0.0,0.0,,,
+probe,mortgage,linear,2.0,2.904149937359671,4.217532586486467,1.4677139768545209,0.0,0.0,,,
+probe,mortgage,linear,2.0,3.148416728152398,4.960202167161096,1.6185830196025777,0.0,0.0,,,
+probe,boston_housing,linear,12.0,1.3520136007298742,0.16837960058355547,0.13424224549151464,0.0,9.0,0.0,1.0,3.862365679999999
+probe,boston_housing,linear,11.0,1.007249276669714,0.09299927052223386,0.09963742976493017,0.0,9.0,,,
+probe,boston_housing,linear,13.0,1.3172484222248864,0.15939193355054326,0.12997258850868726,0.0,9.0,,,
+probe,boston_housing,linear,12.0,1.2995529090837012,0.15490400031856688,0.12785010134920172,0.0,9.0,,,
+probe,boston_housing,linear,12.0,1.0059221718260045,0.09273649552678494,0.09945072011744727,0.0,9.0,,,
 roar,adult,linear,5.0,9.618374680646278,18.55040773192179,2.03277587890625,1.0,5.0,0.05999999999999994,1.0,1.0735000799999999
 roar,adult,linear,5.0,8.90917690170632,15.907446315902012,1.8774079084396362,1.0,5.0,,,
 roar,adult,linear,5.0,14.060748848458747,39.60930977247996,2.9490909576416016,1.0,5.0,,,

experiments/run_experiment.py

@@ -164,6 +164,8 @@ def initialize_recourse_method(
         return Revise(mlmodel, data, hyperparams)
     elif "wachter" in method:
         return Wachter(mlmodel, hyperparams)
+    elif method == "probe":
+        return Probe(mlmodel, hyperparams)
     elif method == "roar":
         return Roar(mlmodel, hyperparams)
     else:
@@ -195,7 +197,7 @@ def create_parser():
     -r, --recourse_method: Specifies recourse methods for the experiment.
         Default: ["dice", "cchvae", "cem", "cem_vae", "clue", "cruds", "face_knn", "face_epsilon", "gs", "mace", "revise", "wachter"].
         Choices: ["dice", "ar", "causal_recourse", "cchvae", "cem", "cem_vae", "claproar", "clue", "cruds", "face_knn", "face_epsilon", "feature_tweak",
-            "focus", "gravitational", "greedy", "gs", "mace", "revise", "wachter", "roar"].
+            "focus", "gravitational", "greedy", "gs", "mace", "revise", "wachter", "roar", "probe"].
     -n, --number_of_samples: Specifies the number of instances per dataset.
         Default: 20.
     -s, --train_split: Specifies the split of the available data used for training.
@@ -284,6 +286,7 @@ def create_parser():
             "mace",
             "revise",
             "wachter",
+            "probe",
             "roar",
         ],
         help="Recourse methods for experiment",
@@ -364,6 +367,7 @@ def create_parser():
         "gravitational",
         "wachter",
         "revise",
+        "probe",
         "roar",
     ]
     sklearn_methods = ["feature_tweak", "focus", "mace"]

methods/__init__.py

@@ -17,6 +17,7 @@
     Gravitational,
     Greedy,
     GrowingSpheres,
+    Probe,
     Revise,
     Roar,
     Wachter,

methods/catalog/__init__.py

@@ -14,6 +14,7 @@
 from .greedy import Greedy
 from .growing_spheres import GrowingSpheres
 from .mace import MACE
+from .probe import Probe
 from .revise import Revise
 from .roar import Roar
 from .wachter import Wachter

methods/catalog/probe/__init__.py

@@ -0,0 +1,4 @@
+# flake8: noqa
+
+from .model import Probe
+from .reproduce import test_probe

methods/catalog/probe/library/__init__.py

@@ -0,0 +1,3 @@
+# flake8: noqa
+
+from .probe import probe_recourse

methods/catalog/probe/library/probe.py

@@ -0,0 +1,286 @@
+import datetime
+from typing import List, Optional
+
+import numpy as np
+import torch
+import torch.distributions.normal as normal_distribution
+import torch.optim as optim
+from torch import nn
+from torch.autograd import Variable
+from torch.distributions.multivariate_normal import MultivariateNormal
+
+from methods.processing import reconstruct_encoding_constraints
+
+"""
+This file contains the implementation of the Probe method, along with required helper functions
+"""
+
+DECISION_THRESHOLD = 0.5
+
+# Mean and variance for rectified normal distribution:
+# see in here : http://journal-sfds.fr/article/view/669
+
+
+def compute_jacobian(inputs, output):
+    """
+    :param inputs: Batch X Size (e.g. Depth X Width X Height)
+    :param output: Batch X Classes
+    :return: jacobian: Batch X Classes X Size
+    """
+    assert inputs.requires_grad
+    grad = gradient(output, inputs)
+    return grad
+
+
+def gradient(y, x, grad_outputs=None):
+    """Compute dy/dx @ grad_outputs"""
+    if grad_outputs is None:
+        grad_outputs = torch.tensor(1, device=y.device)
+    grad = torch.autograd.grad(y, [x], grad_outputs=grad_outputs, create_graph=True)[0]
+    return grad
+
+
+def compute_invalidation_rate_closed(torch_model, x, sigma2):
+    # Compute input into CDF
+    prob = torch_model(x)
+    logit_x = torch.log(prob[0][1] / prob[0][0])
+    Sigma2 = sigma2 * torch.eye(x.shape[0])
+    jacobian_x = compute_jacobian(x, logit_x).reshape(-1)
+    denom = torch.sqrt(sigma2) * torch.norm(jacobian_x, 2)
+    arg = logit_x / denom
+
+    # Evaluate Gaussian cdf
+    normal = normal_distribution.Normal(loc=0.0, scale=1.0)
+    normal_cdf = normal.cdf(arg)
+
+    # Get invalidation rate
+    ir = 1 - normal_cdf
+
+    return ir
+
+
+def perturb_sample(x, n_samples, sigma2):
+    # stack copies of this sample, i.e. n rows of x.
+    X = x.repeat(n_samples, 1)
+    # sample normal distributed values
+    Sigma = torch.eye(x.shape[1]) * sigma2
+    eps = MultivariateNormal(
+        loc=torch.zeros(x.shape[1]), covariance_matrix=Sigma
+    ).sample((n_samples,))
+
+    return X + eps
+
+
+def reparametrization_trick(mu, sigma2, n_samples):
+    # var = torch.eye(mu.shape[1]) * sigma2
+    std = torch.sqrt(sigma2)
+    epsilon = MultivariateNormal(
+        loc=torch.zeros(mu.shape[1]), covariance_matrix=torch.eye(mu.shape[1])
+    )
+    epsilon = epsilon.sample((n_samples,))  # standard Gaussian random noise
+    ones = torch.ones_like(epsilon)
+    random_samples = mu.reshape(-1) * ones + std * epsilon
+
+    return random_samples
+
+
+def compute_invalidation_rate(torch_model, random_samples):
+    yhat = torch_model(random_samples)[:, 1]
+    hat = (yhat > 0.5).float()
+    ir = 1 - torch.mean(hat, 0)
+    return ir
+
+
+def probe_recourse(
+    torch_model,
+    x: np.ndarray,
+    cat_feature_indices: List[int],
+    binary_cat_features: bool = True,
+    feature_costs: Optional[List[float]] = None,
+    lr: float = 0.07,
+    lambda_param: float = 5,
+    y_target: List[int] = [0.45, 0.55],
+    n_iter: int = 500,
+    t_max_min: float = 1.0,
+    norm: int = 1,
+    clamp: bool = False,
+    loss_type: str = "MSE",
+    invalidation_target: float = 0.45,
+    inval_target_eps: float = 0.005,
+    noise_variance: float = 0.01,
+) -> np.ndarray:
+    """
+    Generates counterfactual example according to Wachter et.al for input instance x
+
+    Parameters
+    ----------
+    torch_model: black-box-model to discover
+    x: factual to explain
+    cat_feature_indices: list of positions of categorical features in x
+    binary_cat_features: If true, the encoding of x is done by drop_if_binary
+    feature_costs: List with costs per feature
+    lr: learning rate for gradient descent
+    lambda_param: weight factor for feature_cost
+    y_target: List of one-hot-encoded target class
+    n_iter: maximum number of iteration
+    t_max_min: maximum time of search
+    norm: L-norm to calculate cost
+    clamp: If true, feature values will be clamped to (0, 1)
+    loss_type: String for loss function (MSE or BCE)
+    Invalidation_target: target invalidation rate
+    inval_target_eps: epsilon for invalidation rate
+    noise_variance: variance of the normal distribution for sampling
+
+    Returns
+    -------
+    Counterfactual example as np.ndarray
+    """
+    # device = "cpu"  # for simplicity and to avoid Runtime error.
+    device = "cuda" if torch.cuda.is_available() else "cpu"
+
+    torch_model = torch_model.to(device)
+    # returns counterfactual instance
+    torch.manual_seed(0)
+    noise_variance = torch.tensor(noise_variance)
+
+    # if feature_costs is not None:
+    #     feature_costs = torch.from_numpy(feature_costs).float().to(device)
+
+    # print("x:", x)
+
+    x = torch.from_numpy(x).float().to(device)
+    y_target = torch.tensor(y_target).float().to(device)
+    lamb = torch.tensor(lambda_param).float().to(device)
+    # x_new is used for gradient search in optimizing process
+    x_new = Variable(x.clone(), requires_grad=True)
+    # x_new_enc is a copy of x_new with reconstructed encoding constraints of x_new
+    # such that categorical data is either 0 or 1
+
+    # x_new_enc = reconstruct_encoding_constraints( #TODO: check if this is needed here, i believe that the encoding is done in the model prediction
+    #     x_new, cat_feature_indices, binary_cat_features
+    # )
+
+    optimizer = optim.Adam([x_new], lr, amsgrad=True)
+    softmax = nn.Softmax()
+
+    if loss_type == "MSE":
+        loss_fn = torch.nn.MSELoss()
+        f_x_new = softmax(torch_model(x_new))[:, 1]
+    else:
+        loss_fn = torch.nn.BCELoss()
+        f_x_new = torch_model(x_new)[:, 1]
+
+    t0 = datetime.datetime.now()
+    t_max = datetime.timedelta(minutes=t_max_min)
+
+    costs = []
+    ces = []
+
+    random_samples = reparametrization_trick(x_new, noise_variance, n_samples=1000)
+    invalidation_rate = compute_invalidation_rate(torch_model, random_samples)
+
+    while (f_x_new <= DECISION_THRESHOLD) or (
+        invalidation_rate > invalidation_target + inval_target_eps
+    ):
+        # it = 0
+        for it in range(n_iter):
+            # while invalidation_target >= 0.5 and it < n_iter:
+
+            optimizer.zero_grad()
+            # x_new_enc = reconstruct_encoding_constraints(
+            #     x_new, cat_feature_indices, binary_cat_features
+            # )
+            # use x_new_enc for prediction results to ensure constraints
+            # f_x_new = softmax(torch_model(x_new))[:, 1]
+            f_x_new_binary = torch_model(x_new).squeeze(axis=0)
+
+            cost = (
+                torch.dist(x_new, x, norm)
+                # if feature_costs is None
+                # else torch.norm(feature_costs * (x_new - x), norm)
+            )
+
+            # Compute Invalidation loss
+            # output_mean, output_std = compute_output_dist_suff_statistics(torch_model, x_new,
+            #                                                              noise_variance=noise_variance)
+
+            # normal = normal_distribution.Normal(loc=0.0, scale=1.0)
+            # ratio = torch.divide(output_mean, output_std)
+            # normal_cdf = normal.cdf(ratio)
+            # invalidation_rate = 1 - normal_cdf
+
+            # invalidation_rate = compute_invalidation_rate(torch_model, random_samples)
+            invalidation_rate_c = compute_invalidation_rate_closed(
+                torch_model, x_new, noise_variance
+            )
+
+            # Compute & update losses
+            loss_invalidation = invalidation_rate_c - invalidation_target
+            # Hinge loss
+            loss_invalidation[loss_invalidation < 0] = 0
+
+            loss = (
+                3 * loss_invalidation + loss_fn(f_x_new_binary, y_target) + lamb * cost
+            )
+            loss.backward()
+            optimizer.step()
+
+            random_samples = reparametrization_trick(
+                x_new, noise_variance, n_samples=10000
+            )
+            invalidation_rate = compute_invalidation_rate(torch_model, random_samples)
+
+            # x_pertub = perturb_sample(x_new, sigma2=noise_variance, n_samples=10000)
+            # pred = 1 - torch_model(x_pertub)[:, 1]
+            # invalidation_rate_empirical = torch.mean(pred)
+
+            # print('-----------------------------------------')
+            # print('IR empirical', invalidation_rate_empirical)
+            # print('IR from loss', invalidation_rate)
+            # print('IR loss', loss_invalidation)
+
+            # clamp potential CF
+            if clamp:
+                x_new.clone().clamp_(0, 1)
+            # it += 1
+
+            # x_new_enc = reconstruct_encoding_constraints(
+            #     x_new, cat_feature_indices, binary_cat_features
+            # )
+            # f_x_new = torch_model(x_new_enc)[:, 1]
+            f_x_new = torch_model(x_new)[:, 1]
+
+        if (f_x_new > DECISION_THRESHOLD) and (
+            invalidation_rate < invalidation_target + inval_target_eps
+        ):
+            print("--------------------------------------")
+            print("invalidation rate:", invalidation_rate)
+            # print('emp invalidation rate', invalidation_rate_empirical)
+            print("cost:", cost)
+            print("classifier output:", f_x_new_binary)
+
+            costs.append(cost)
+            ces.append(x_new)
+
+            break
+
+        lamb -= 0.10
+
+        if datetime.datetime.now() - t0 > t_max:
+            print("Timeout")
+            break
+
+    if not ces:
+        print(
+            "No Counterfactual Explanation Found at that Target Rate - Try Different Target"
+        )
+        return x_new.cpu().detach().numpy().squeeze(axis=0)
+    else:
+        print("Counterfactual Explanation Found")
+        costs = torch.tensor(costs)
+        min_idx = int(torch.argmin(costs).numpy())
+        x_new_enc = ces[min_idx]
+
+    # print("x_prime ", x_new_enc.cpu().detach().numpy().squeeze(axis=0))
+
+    return x_new_enc.cpu().detach().numpy().squeeze(axis=0)