experiment1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from openxai.model import LoadModel\n",
    "from openxai.dataloader import return_loaders\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.genetic import GeneticAlgorithm\n",
    "from src.explainer import Explainer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_name = 'pima' # gaussian, heloc, pima\n",
    "model_name = 'ann'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loader_train, loader_test = return_loaders(data_name=data_name, download=True, batch_size=128)\n",
    "# X = loader_test.dataset.data # gaussian, heloc\n",
    "X = loader_train.dataset.data # pima as test is too small to estimate cd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = LoadModel(data_name=data_name, ml_model=model_name, pretrained=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def partial_dependence(f, xs, X, s):\n",
    "    X_raw = X.copy()\n",
    "    X_raw[:, s] = xs\n",
    "    return f.predict_proba(X_raw)[:, 1].mean()\n",
    "\n",
    "def conditional_dependence(f, xs, X, s):\n",
    "    X_raw = X.copy()\n",
    "    X_s = X_raw[:, s]\n",
    "    epsilon = X_s.ptp() / 18\n",
    "    X_cond = X_raw[(X_s > xs - epsilon) & (X_s < xs + epsilon), :]\n",
    "    if X_cond.shape[0] == 0:\n",
    "        return partial_dependence(f, xs, X, s)\n",
    "    else:\n",
    "        X_cond[:, s] = xs\n",
    "        return f.predict_proba(X_cond)[:, 1].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "X_s = X[:, i]\n",
    "print(partial_dependence(model, np.quantile(X_s, 0.25), X, i))\n",
    "print(conditional_dependence(model, np.quantile(X_s, 0.25), X, i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_dtv_nd(X, Y, method=\"histogram\", nbins=None):\n",
    "    if method == \"histogram\":\n",
    "        if nbins is None:\n",
    "            nbins = int(np.power(X.shape[0] / 10, 1 / X.shape[1]))\n",
    "            if nbins < 5:\n",
    "                nbins = 5\n",
    "        XY = np.concatenate((X, Y))\n",
    "        _, edges = np.histogramdd(XY, density=True, bins=nbins)\n",
    "        X_density = np.histogramdd(X, density=False, bins=edges)[0] / X.shape[0]\n",
    "        Y_density = np.histogramdd(Y, density=False, bins=edges)[0] / Y.shape[0]\n",
    "    elif method == \"kernel\":\n",
    "        X_kernel = stats.gaussian_kde(X.T)\n",
    "        Y_kernel = stats.gaussian_kde(Y.T)\n",
    "        XY = np.concatenate((X.T, Y.T), axis=1)\n",
    "        X_density = X_kernel.pdf(XY) / X_kernel.pdf(XY).sum()\n",
    "        Y_density = Y_kernel.pdf(XY) / Y_kernel.pdf(XY).sum()\n",
    "    return np.sum(np.abs(X_density - Y_density)) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = {}\n",
    "for s in range(X.shape[1]):\n",
    "    temp[s] = []\n",
    "    X_s = X[:, i]\n",
    "    for xs in np.linspace(X_s.min(), X_s.max(), 18):\n",
    "        temp[s] += [partial_dependence(model, xs, X, s)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_pd = pd.DataFrame(temp)\n",
    "feature_importance_rank = temp_pd.var(axis=0).argsort().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = feature_importance_rank[len(feature_importance_rank)//2]\n",
    "xs = X[:, s].mean()\n",
    "corr_X = np.corrcoef(X.T)\n",
    "np.fill_diagonal(corr_X, 0)\n",
    "K = 4\n",
    "# features_to_change = corr_X[:, s].argsort()[-K:]\n",
    "features_to_change = feature_importance_rank[-np.array(range(1, K+1))]\n",
    "exp = Explainer(model, xs, X, s, \"pd\")\n",
    "alg = GeneticAlgorithm(exp, np.setdiff1d(list(range(X.shape[1])), features_to_change), stop_iter=100)\n",
    "alg.attack(max_iter=100, random_state=123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "\n",
    "bounds = pd.DataFrame({'s': [], 'xs': [], 'dfe': [], 'label': []})\n",
    "\n",
    "for s_name in ['most important', 'median important', 'least important']:\n",
    "    features_to_change = feature_importance_rank[-np.array(range(1, K+1))]\n",
    "    if s_name == \"most important\":\n",
    "        s = feature_importance_rank[-1]\n",
    "        features_to_change = feature_importance_rank[-np.array(range(2, K+2))]\n",
    "    elif s_name == \"median important\":\n",
    "        s = feature_importance_rank[len(feature_importance_rank)//2]\n",
    "    elif s_name == \"least important\":\n",
    "        s = feature_importance_rank[0]\n",
    "    X_s = X[:, s]                        \n",
    "\n",
    "    for q in [0.2, 0.5, 0.8]:\n",
    "        xs = np.quantile(X_s, q)\n",
    "        pd1 = partial_dependence(model, xs, X, s)\n",
    "        cd1 = conditional_dependence(model, xs, X, s)\n",
    "        X2 = X.copy()\n",
    "        ###\n",
    "        temp = X.copy()\n",
    "        temp[:, s] = xs\n",
    "        ###\n",
    "        B_xs = model.predict_proba(temp)[:, 1].max()\n",
    "        print(f'----- s = {s_name} | x_s = q_{q} | B_xs = {B_xs} | pd_s = {pd1} | cd_s = {cd1} -----')\n",
    "\n",
    "        bounds = pd.concat([bounds, pd.DataFrame({\n",
    "            's': [s_name]*2, \n",
    "            'xs': ['q_'+str(q)]*2, \n",
    "            'dfe': [np.max([np.abs(B_xs - pd1), pd1]), np.max([np.abs(B_xs - cd1), cd1])], \n",
    "            'label': ['marginal', 'conditional']\n",
    "        })])\n",
    "bounds = bounds.assign(xs=bounds[\"xs\"].str.replace(\"_\", \"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds.to_csv(f'results/{model_name}_{data_name}_remark2.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "B = model.predict_proba(X)[:, 1].max()\n",
    "\n",
    "result = pd.DataFrame({'s': [], 'xs': [], 'dtv': [], 'dfe': [], 'label': []})\n",
    "\n",
    "ITER = 200\n",
    "STOP_ITER = 30\n",
    "POP_COUNT = 100\n",
    "TOP_SURVIVORS = 5\n",
    "\n",
    "corr_X = np.corrcoef(X.T)\n",
    "np.fill_diagonal(corr_X, 0)\n",
    "K = 2\n",
    "\n",
    "for s_name in ['most important', 'median important', 'least important']:\n",
    "    features_to_change = feature_importance_rank[-np.array(range(1, K+1))]\n",
    "    if s_name == \"most important\":\n",
    "        s = feature_importance_rank[-1]\n",
    "        features_to_change = feature_importance_rank[-np.array(range(2, K+2))]\n",
    "    elif s_name == \"median important\":\n",
    "        s = feature_importance_rank[len(feature_importance_rank)//2]\n",
    "    elif s_name == \"least important\":\n",
    "        s = feature_importance_rank[0]\n",
    "    X_s = X[:, s]                        \n",
    "    epsilon = X_s.ptp() / 18\n",
    "    # features_to_change = corr_X[:, s].argsort()[-K:]\n",
    "    for q in [0.2, 0.5, 0.8]:\n",
    "        xs = np.quantile(X_s, q)\n",
    "        pd1 = partial_dependence(model, xs, X, s)\n",
    "        cd1 = conditional_dependence(model, xs, X, s)\n",
    "        X2 = X.copy()\n",
    "        ###\n",
    "        temp = X.copy()\n",
    "        temp[:, s] = xs\n",
    "        ###\n",
    "        print(f'----- s = {s_name} | x_s = q_{q} | B_xs = {model.predict_proba(temp)[:, 1].max()} -----')\n",
    "        for perturbation in ['random', 'adv']:\n",
    "            if perturbation == \"adv\":\n",
    "                for i, sigma in enumerate([0.01, 0.05, 0.1, 0.25, 0.33]): \n",
    "                    TARGET_PD = 0\n",
    "                    TARGET_CD = 0\n",
    "                    for seed in range(5):\n",
    "                        exp_pd = Explainer(model, xs, X, s, \"pd\")\n",
    "                        alg_pd = GeneticAlgorithm(exp_pd, np.setdiff1d(list(range(X.shape[1])), features_to_change), \n",
    "                                                  std_ratio=sigma, pop_count=POP_COUNT, top_survivor=TOP_SURVIVORS, stop_iter=STOP_ITER)\n",
    "                        alg_pd.attack(max_iter=ITER, random_state=seed, cache_data=True, target=TARGET_PD)\n",
    "                        X2 = alg_pd.get_best_data()\n",
    "                        pd2 = partial_dependence(model, xs, X2, s)\n",
    "                        d_pd = np.abs(pd1 - pd2)\n",
    "                        dtv_pd = distance_dtv_nd(X[:, features_to_change], X2[:, features_to_change])\n",
    "\n",
    "                        exp_cd = Explainer(model, xs, X, s, \"cd\")\n",
    "                        alg_cd = GeneticAlgorithm(exp_cd, np.setdiff1d(list(range(X.shape[1])), features_to_change), \n",
    "                                                  std_ratio=sigma, pop_count=POP_COUNT, top_survivor=TOP_SURVIVORS, stop_iter=STOP_ITER)\n",
    "                        alg_cd.attack(max_iter=ITER, random_state=seed, cache_data=True, target=TARGET_CD)\n",
    "                        X2 = alg_cd.get_best_data()\n",
    "                        cd2 = conditional_dependence(model, xs, X2, s)\n",
    "                        d_cd = np.abs(cd1 - cd2)\n",
    "                        X_cond = X[(X_s > xs - epsilon) & (X_s < xs + epsilon), :]\n",
    "                        X2_cond = X2[(X2[:, s] > xs - epsilon) & (X2[:, s] < xs + epsilon), :]\n",
    "                        dtv_cd = distance_dtv_nd(X_cond[:, features_to_change], X2_cond[:, features_to_change])\n",
    "\n",
    "                        if seed == 0: # check if alternative target is better\n",
    "                            alg_pd = GeneticAlgorithm(exp_pd, np.setdiff1d(list(range(X.shape[1])), features_to_change), \n",
    "                                                      std_ratio=sigma, pop_count=POP_COUNT, top_survivor=TOP_SURVIVORS, stop_iter=STOP_ITER)\n",
    "                            alg_pd.attack(max_iter=ITER, random_state=seed, cache_data=True, target=1 - TARGET_PD)\n",
    "                            X2 = alg_pd.get_best_data()\n",
    "                            pd2 = partial_dependence(model, xs, X2, s)\n",
    "                            d_pd2 = np.abs(pd1 - pd2) \n",
    "                            if d_pd < d_pd2: # if better: update \n",
    "                                d_pd = d_pd2\n",
    "                                dtv_pd = distance_dtv_nd(X[:, features_to_change], X2[:, features_to_change])\n",
    "                                TARGET_PD = 1 - TARGET_PD\n",
    "\n",
    "                            alg_cd = GeneticAlgorithm(exp_cd, np.setdiff1d(list(range(X.shape[1])), features_to_change), \n",
    "                                                      std_ratio=sigma, pop_count=POP_COUNT, top_survivor=TOP_SURVIVORS, stop_iter=STOP_ITER)\n",
    "                            alg_cd.attack(max_iter=ITER, random_state=seed, cache_data=True, target=1 - TARGET_CD)\n",
    "                            X2 = alg_cd.get_best_data()\n",
    "                            cd2 = conditional_dependence(model, xs, X2, s)\n",
    "                            d_cd2 = np.abs(cd1 - cd2)\n",
    "                            if d_cd < d_cd2: # if better: update \n",
    "                                d_cd = d_cd2\n",
    "                                X_cond = X[(X_s > xs - epsilon) & (X_s < xs + epsilon), :]\n",
    "                                X2_cond = X2[(X2[:, s] > xs - epsilon) & (X2[:, s] < xs + epsilon), :]\n",
    "                                dtv_cd = distance_dtv_nd(X_cond[:, features_to_change], X2_cond[:, features_to_change])\n",
    "                                TARGET_CD = 1 - TARGET_CD\n",
    "\n",
    "                        result = pd.concat([result, pd.DataFrame({\n",
    "                            's': [s_name]*2, \n",
    "                            'xs': ['q_'+str(q)]*2, \n",
    "                            'dtv': [dtv_pd, dtv_cd],\n",
    "                            'dfe': [d_pd, d_cd], \n",
    "                            'label': ['adversarial:marginal', 'adversarial:conditional']\n",
    "                        })])\n",
    "\n",
    "            if perturbation == \"random\":\n",
    "                for i, sigma in enumerate([0.01, 0.05, 0.1, 0.12, 0.25]): \n",
    "                    for seed in range(5):\n",
    "                        np.random.seed(seed)\n",
    "                        X2[:, features_to_change] += np.random.normal(0, sigma, size=(X2.shape[0], K))\n",
    "                        pd2 = partial_dependence(model, xs, X2, s)\n",
    "                        cd2 = conditional_dependence(model, xs, X2, s)\n",
    "                        dtv_pd = distance_dtv_nd(X[:, features_to_change], X2[:, features_to_change])\n",
    "\n",
    "                        X_cond = X[(X_s > xs - epsilon) & (X_s < xs + epsilon), :]\n",
    "                        X2_cond = X2[(X2[:, s] > xs - epsilon) & (X2[:, s] < xs + epsilon), :]\n",
    "                        dtv_cd = distance_dtv_nd(X_cond[:, features_to_change], X2_cond[:, features_to_change])\n",
    "\n",
    "                        result = pd.concat([result, pd.DataFrame({\n",
    "                            's': [s_name]*4, \n",
    "                            'xs': ['q_'+str(q)]*4, \n",
    "                            'dtv': [dtv_pd, dtv_cd, dtv_pd, dtv_cd],\n",
    "                            'dfe': [2*B*dtv_pd, 2*B*dtv_cd, np.abs(pd1 - pd2), np.abs(cd1 - cd2)], \n",
    "                            'label': ['theoretical:marginal', 'theoretical:conditional', 'random:marginal', 'random:conditional']\n",
    "                        })])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = result.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv(f'results/{model_name}_{data_name}_k2.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.assign(xs=df[\"xs\"].str.replace(\"_\", \"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.relplot(\n",
    "    data=df.loc[(df.dfe < B+.01) & (df.dtv < 0.7), :], \n",
    "    x=\"dtv\", y=\"dfe\",\n",
    "    col=\"xs\", row=\"s\", \n",
    "    hue=\"label\", style=\"label\",\n",
    "    kind=\"line\", linewidth=2,\n",
    "    height=1.8,\n",
    "    aspect=1.1,\n",
    "    palette=[\"#4285F4\", \"#DB4437\", \"#4285F4\", \"#DB4437\", \"#4285F4\", \"#DB4437\"], # [\"black\", \"#0F9D58\", \"#F4B400\", \"#4285F4\", \"#DB4437\"]\n",
    "    dashes=[(1, 1), (1, 1), (2, 2), (2, 2), \"\", \"\"]\n",
    ")\n",
    "g._legend.set_title(\"Input perturbation\")\n",
    "g.set_axis_labels(\"$d_{\\mathrm{TV}}(p_{\\mathbf{X}}, p'_{\\mathbf{X}})$\", \n",
    "                  \"$| g_s(\\mathbf{x}_s; f, p_{\\mathbf{X}}) - g_s(\\mathbf{x}_s; f, p'_{\\mathbf{X}}) |$\")\n",
    "for i, ax in enumerate(g.axes):\n",
    "    if i == 1:\n",
    "        ax[0].set_ylabel(\"\")\n",
    "g.fig.subplots_adjust(top=0.9)\n",
    "sns.move_legend(g, \"upper left\", bbox_to_anchor=(.15, 0.65), frameon=True, ncol=3) # (.45, .38)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"../figures/exp1_heloc.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = pd.read_csv(f'results/{model_name}_{data_name}_k2.csv') # gaussian, heloc, pima\n",
    "df2 = df2.assign(xs=df2[\"xs\"].str.replace(\"_\", \"\"))\n",
    "\n",
    "bounds = pd.read_csv(f'results/{model_name}_{data_name}_remark2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, row in bounds.iterrows():\n",
    "    ix = (df2['s'] == row['s']) &\\\n",
    "        (df2['xs'] == row['xs']) &\\\n",
    "        (df2['label'].str.split(\":\", expand=True).iloc[:, 0] == \"theoretical\") &\\\n",
    "        (df2['label'].str.split(\":\", expand=True).iloc[:, 1] == row['label']) &\\\n",
    "        (df2['dfe'] > row['dfe'])\n",
    "    df2.loc[ix, \"dfe\"] = row['dfe']\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.relplot(\n",
    "    data=df2,\n",
    "    x=\"dtv\", y=\"dfe\",\n",
    "    col=\"xs\", row=\"s\", \n",
    "    hue=\"label\", style=\"label\",\n",
    "    kind=\"line\", linewidth=2,\n",
    "    height=1.8,\n",
    "    aspect=1.1,\n",
    "    palette=[\"#4285F4\", \"#DB4437\", \"#4285F4\", \"#DB4437\", \"#4285F4\", \"#DB4437\"], # [\"black\", \"#0F9D58\", \"#F4B400\", \"#4285F4\", \"#DB4437\"]\n",
    "    dashes=[(1, 1), (1, 1), (2, 2), (2, 2), \"\", \"\"]\n",
    ")\n",
    "g._legend.set_title(\"Input perturbation\")\n",
    "g.set_axis_labels(\"$d_{\\mathrm{TV}}(p_{\\mathbf{X}}, p'_{\\mathbf{X}})$\", \n",
    "                  \"$| g_s(\\mathbf{x}_s; f, p_{\\mathbf{X}}) - g_s(\\mathbf{x}_s; f, p'_{\\mathbf{X}}) |$\")\n",
    "for i, ax in enumerate(g.axes):\n",
    "    if i == 1:\n",
    "        ax[0].set_ylabel(\"\")\n",
    "# g.fig.subplots_adjust(top=1.5)\n",
    "sns.move_legend(g, \"upper left\", bbox_to_anchor=(.15, 1.13), frameon=True, ncol=3) # (.45, .38)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"../figures/exp1_pima_tighter.pdf\", bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "robustfe",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}