basic_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "import pathlib\n",
    "from os.path import join\n",
    "path_to_file = str(pathlib.Path().resolve())\n",
    "dir_path = join(path_to_file)\n",
    "\n",
    "sys.path.append(join(dir_path, \"HelperFiles\"))\n",
    "import helper\n",
    "import retrospective\n",
    "import top_k \n",
    "import train_models\n",
    "import load_data \n",
    "\n",
    "from os.path import join\n",
    "data_path = join(dir_path, \"Experiments\", \"Data\")\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from scipy.stats import binom\n",
    "import importlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dataset = \"census\" #breast_cancer\n",
    "# X_train, y_train, X_test, y_test, mapping_dict = load_data.load_data(data_path, dataset)\n",
    "# model = train_models.train_model(X_train, y_train, \"nn\", lime=False)\n",
    "\n",
    "# print(\"Class imbalance: {}%\".format(round(100*(max(np.mean(y_test), 1-np.mean(y_test))))))\n",
    "# Y_preds = (model(X_test) > 0.5).astype(\"int\")\n",
    "# print(\"NN {}% accuracy\".format(round(np.mean(Y_preds == y_test)*100)))\n",
    "# d = len(mapping_dict) if mapping_dict is not None else X_train.shape[1]\n",
    "# print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ongoing runs\n",
    "## Retrospective rank verification\n",
    "### Shapley Sampling and KernelSHAP, various alphas. Both always control FWER!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "results_path = join(dir_path, \"Experiments\", \"Results\")\n",
    "retro_path = join(results_path, \"Retrospective\")\n",
    "method = \"ss\"\n",
    "dataset = \"census\"\n",
    "with open(join(retro_path, method+\"_\"+dataset), 'rb') as f:\n",
    "    retro_results = pickle.load(f)\n",
    "shap_vals = retro_results[\"shap_vals\"]\n",
    "N_verified_all = retro_results[\"N_verified\"]\n",
    "shap_vars = retro_results[\"shap_vars\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alphas:\t [0.05, 0.1, 0.2]\n",
      "(30, 50, 3)\n",
      "Max:\t [0.02 0.02 0.1 ]\n",
      "Mean:\t [0.004 0.008 0.021]\n"
     ]
    }
   ],
   "source": [
    "# Number of stable ranks from retrospective procedure. \n",
    "importlib.reload(helper)\n",
    "alphas = [0.05, 0.1, 0.2]\n",
    "print(\"Alphas:\\t\", alphas)\n",
    "all_ranks = helper.shap_vals_to_ranks(shap_vals, abs=True)\n",
    "\n",
    "N_pts, N_runs, N_alphas = N_verified_all.shape\n",
    "print(N_verified_all.shape)\n",
    "avg_shap = np.mean(shap_vals, axis=1)\n",
    "avg_ranks = np.array([helper.get_ranking(avg_shap[i], abs=True) for i in range(N_pts)])\n",
    "\n",
    "fwers = helper.calc_all_retro_fwers(N_verified_all, all_ranks, avg_ranks)\n",
    "print(\"Max:\\t\", np.round(np.nanmax(fwers, axis=1), 3))\n",
    "print(\"Mean:\\t\", np.round(np.nanmean(fwers, axis=1), 3))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAADcCAYAAABtaehmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAf5UlEQVR4nO3de1hU1f4G8HcAuemAchMvKAgYiHgJJA8aEFpGJpKImlFGmfYcDcnLibwEnlAgM+2gpqiheamIkJTykCkGlh01o8I0UMF4NAE7wnDxoDL794c/psYBBQS2s3g/z8Mfs/aevb57M8zLXmvvGYUkSRKIiIhIGAZyF0BERERti+FOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U7COnz4MBQKBdLS0uQupVlKS0sxefJkWFtbQ6FQYO3atXKXBABQKBSIjY3Vajt+/Dh8fX3RtWtXKBQK5OXlITY2FgqFok37DggIQEBAQJtu81417OeVK1fkLoWoSUZyF0D6bdu2bYiIiICJiQnOnTuHPn36aC0PCAjAlStXkJ+fL1OF+uPVV19FVlYWYmJiYG9vD29vb7lLatSNGzcQFhYGU1NTrFmzBubm5ujfv7/cZRHRXzDcqU3U1dUhISEBSUlJcpeitw4dOoSJEydi4cKFcpei5dq1azAy+vOt4ty5c7hw4QI2b96MmTNnatqXLl2K6OhoOUokottwWJ7axLBhw7B582ZcunRJ7lI6XE1NTZtsp6ysDN27d2+Tbd0rtVqN//3vfwAAU1NTrXAvKysDAJ1ajYyMYGpq2mE13slf6yfqjBju1CYWL16M+vp6JCQk3HG94uJiKBQKbNu2TWfZ7XO7DXObBQUFCA8Ph6WlJWxtbbFs2TJIkoSSkhJMnDgRFhYWsLe3x+rVqxvts76+HosXL4a9vT26du2K4OBglJSU6Kz3n//8B48//jgsLS1hbm4Of39/fPPNN1rrNNT0yy+/YPr06ejRowdGjx59x30+f/48wsLCYGVlBXNzc4wcORKff/65Zvm2bdugUCggSRLWr18PhULR5Nz1jRs3YGVlhYiICJ1lKpUKpqamWmf+dXV1iImJgYuLC0xMTODg4IB//OMfqKur03quQqHA3LlzsWvXLnh4eMDExAT//ve/Ncsafi/PP/88/P39AQBhYWFQKBSaOfGm5tx37twJLy8vmJmZwcrKCtOmTWv0+CcnJ8PZ2RlmZmbw8fFBbm7uHY6qtjvV//bbb8PX1xfW1tYwMzODl5dXo9dhNGwjIyMDgwcPhomJCTw8PDTbuZMLFy7AxcUFgwcPRmlpKQCgsLAQoaGhsLe3h6mpKfr27Ytp06ahsrKy2fvV4MaNG1i+fDlcXV1hamoKa2trjB49GgcOHNCsc/nyZURERKBv374wMTFBr169MHHiRBQXF7e4P9J/HJanNuHk5ITnnnsOmzdvRnR0NHr37t1m2546dSrc3d2RkJCAzz//HHFxcbCyssKmTZsQGBiIxMRE7Nq1CwsXLsSIESPg5+en9fwVK1ZAoVDgtddeQ1lZGdauXYuxY8ciLy8PZmZmAG4NiQcFBcHLywsxMTEwMDBASkoKAgMDkZubCx8fH61thoWFwdXVFStXrsSdvjW5tLQUvr6+qK2tRWRkJKytrbF9+3YEBwcjLS0NTz31FPz8/LBjxw48++yzePTRR/Hcc881ub0uXbrgqaeeQnp6OjZt2gRjY2PNsoyMDNTV1WHatGkAbp29BgcH48iRI5g1axbc3d3x888/Y82aNSgoKEBGRobWtg8dOoTU1FTMnTsXNjY2cHR01Ol/9uzZ6NOnD1auXInIyEiMGDECPXv2bLLeFStWYNmyZZgyZQpmzpyJ8vJyJCUlwc/PDz/88IPm7H/r1q2YPXs2fH19ERUVhfPnzyM4OBhWVlZwcHBocvvNqf/dd99FcHAwnnnmGVy/fh0fffQRwsLCkJmZifHjx2tt48iRI0hPT8ff//53KJVK/Otf/0JoaCh+++03WFtbN9rvuXPnEBgYCCsrKxw4cAA2Nja4fv06xo0bh7q6Orzyyiuwt7fHxYsXkZmZiYqKClhaWjZrnxrExsYiPj4eM2fOhI+PD1QqFU6cOIGTJ0/i0UcfBQCEhobi1KlTeOWVV+Do6IiysjIcOHAAv/32W6O/SxKcRHQPUlJSJADS8ePHpXPnzklGRkZSZGSkZrm/v7/k4eGheVxUVCQBkFJSUnS2BUCKiYnRPI6JiZEASLNmzdK03bx5U+rbt6+kUCikhIQETfvVq1clMzMzacaMGZq27OxsCYDUp08fSaVSadpTU1MlANK7774rSZIkqdVqydXVVRo3bpykVqs169XW1kpOTk7So48+qlPT008/3azjExUVJQGQcnNzNW1VVVWSk5OT5OjoKNXX12vt/5w5c+66zaysLAmAtG/fPq32J554QhowYIDm8Y4dOyQDAwOtviVJkjZu3CgBkL755hutvg0MDKRTp07p9Hf776XhuH7yySda6zUcmwbFxcWSoaGhtGLFCq31fv75Z8nIyEjTfv36dcnOzk4aNmyYVFdXp1kvOTlZAiD5+/vf5Yjcuf7a2lqtx9evX5cGDx4sBQYG6mzD2NhYOnv2rKbtxx9/lABISUlJOvtZXl4unT59Wurdu7c0YsQI6b///a9mnR9++KHRY9RaQ4cOlcaPH9/k8qtXr0oApFWrVrVJf6T/OCxPbWbAgAF49tlnkZycjN9//73NtvvXi7YMDQ3h7e0NSZLw4osvatq7d++OBx54AOfPn9d5/nPPPQelUql5PHnyZPTq1QtffPEFACAvLw+FhYWYPn06/vjjD1y5cgVXrlxBTU0NxowZg5ycHKjVaq1tvvzyy82q/YsvvoCPj4/W0H23bt0wa9YsFBcX45dffmneQfiLwMBA2NjY4OOPP9a0Xb16FQcOHMDUqVM1bZ988gnc3d3h5uam2acrV64gMDAQAJCdna21XX9/fwwaNKjF9TQlPT0darUaU6ZM0erf3t4erq6umv5PnDiBsrIyvPzyy1ojEc8//3yLznCbqr9hdAa4dZwqKyvx8MMP4+TJkzrrjh07Fs7OzprHQ4YMgYWFRaOvq/z8fPj7+8PR0RFfffUVevTooVnWUHdWVhZqa2ubvQ9N6d69O06dOoXCwsJGl5uZmcHY2BiHDx/G1atX77k/0n8Md2pTS5cuxc2bN+86994S/fr103psaWkJU1NT2NjY6LQ39sbm6uqq9VihUMDFxUUzF9nwhjljxgzY2tpq/WzZsgV1dXU686ROTk7Nqv3ChQt44IEHdNrd3d01y1vKyMgIoaGh+OyzzzRz5+np6bhx44ZWuBcWFuLUqVM6+zRw4EAAf14Y19J9aq7CwkJIkgRXV1edGk6fPq3pv+EY3P576tKlCwYMGNDs/pqqPzMzEyNHjoSpqSmsrKxga2uL9957r9G579tfawDQo0ePRl9XEyZMgFKpRFZWFiwsLHRqmT9/PrZs2QIbGxuMGzcO69evb9V8OwD885//REVFBQYOHAhPT08sWrQIP/30k2a5iYkJEhMTsX//fvTs2RN+fn546623cPny5Vb1R/qPc+7UpgYMGIDw8HAkJyc3eltUUxeK1dfXN7lNQ0PDZrUBuOP8d1MazspXrVqFYcOGNbpOt27dtB7/9WxQDtOmTcOmTZuwf/9+hISEIDU1FW5ubhg6dKhmHbVaDU9PT7zzzjuNbuP2uey23ie1Wg2FQoH9+/c3+vu6/Zjeq8bqz83NRXBwMPz8/LBhwwb06tULXbp0QUpKCnbv3q2zfkteV6Ghodi+fTt27dqF2bNn6yxfvXo1nn/+eXz22Wf48ssvERkZifj4eHz33Xfo27dvi/bNz88P586d02xry5YtWLNmDTZu3KgZ2YqKisKECROQkZGBrKwsLFu2DPHx8Th06BCGDx/eov5I/zHcqc0tXboUO3fuRGJios6yhqHLiooKrfbWnME21+1DmZIk4ezZsxgyZAgAaIZhLSwsMHbs2Dbtu3///vj111912s+cOaNZ3hp+fn7o1asXPv74Y4wePRqHDh3CkiVLtNZxdnbGjz/+iDFjxrT5J8c1h7OzMyRJgpOTk2a0oDENx6CwsFAzZQDcukK8qKhI6x+Wlvr0009hamqKrKwsmJiYaNpTUlJavc0Gq1atgpGRkebiu+nTp+us4+npCU9PTyxduhTffvstRo0ahY0bNyIuLq7F/TXcJREREYHq6mr4+fkhNjZWa9rK2dkZCxYswIIFC1BYWIhhw4Zh9erV2Llz5z3tK+kfDstTm3N2dkZ4eDg2bdqkMyxoYWEBGxsb5OTkaLVv2LCh3er54IMPUFVVpXmclpaG33//HUFBQQAALy8vODs74+2330Z1dbXO88vLy1vd9xNPPIFjx47h6NGjmraamhokJyfD0dGx1XPcBgYGmDx5Mvbt24cdO3bg5s2bWkPyADBlyhRcvHgRmzdv1nn+tWvX2uz+/KZMmjQJhoaGWL58uc6ZryRJ+OOPPwAA3t7esLW1xcaNG3H9+nXNOtu2bdP5J7ClDA0NoVAotEaGiouLde4UaA2FQoHk5GRMnjwZM2bMwN69ezXLVCoVbt68qbW+p6cnDAwMdG5DbI6GY9WgW7ducHFx0WyrtrZW575+Z2dnKJXKVvVH+o9n7tQulixZgh07duDXX3+Fh4eH1rKZM2ciISEBM2fOhLe3N3JyclBQUNButVhZWWH06NGIiIhAaWkp1q5dCxcXF7z00ksAbgXlli1bEBQUBA8PD0RERKBPnz64ePEisrOzYWFhgX379rWq7+joaHz44YcICgpCZGQkrKyssH37dhQVFeHTTz+FgUHr/7+eOnUqkpKSEBMTA09PT808foNnn30WqampePnll5GdnY1Ro0ahvr4eZ86cQWpqKrKystr1I26dnZ0RFxeH119/HcXFxQgJCYFSqURRURH27NmDWbNmYeHChejSpQvi4uIwe/ZsBAYGYurUqSgqKkJKSkqL5twbM378eLzzzjt4/PHHMX36dJSVlWH9+vVwcXHRmrNuLQMDA+zcuRMhISGYMmUKvvjiCwQGBuLQoUOYO3cuwsLCMHDgQNy8eRM7duyAoaEhQkNDNc+PjY3F8uXLkZ2dfcfP0B80aBACAgLg5eUFKysrnDhxAmlpaZg7dy4AoKCgAGPGjMGUKVMwaNAgGBkZYc+ePSgtLdXcGkmdC8Od2oWLiwvCw8Oxfft2nWVvvPEGysvLkZaWhtTUVAQFBWH//v2ws7Nrl1oWL16Mn376CfHx8aiqqsKYMWOwYcMGmJuba9YJCAjA0aNH8eabb2LdunWorq6Gvb09HnrooUbnU5urZ8+e+Pbbb/Haa68hKSkJ//vf/zBkyBDs27dP5x7rlvL19YWDgwNKSkp0ztqBW8GTkZGBNWvW4IMPPsCePXtgbm6OAQMGYN68eXccKm8r0dHRGDhwINasWYPly5cDuDXX/9hjjyE4OFiz3qxZs1BfX49Vq1Zh0aJF8PT0xN69e7Fs2bJ76j8wMBBbt25FQkICoqKi4OTkhMTERBQXF7dJuAO3LvxLS0tDUFAQJk6ciK+++gpDhw7FuHHjsG/fPly8eBHm5uYYOnQo9u/fj5EjR2qeW11dDYVCAXt7+zv2ERkZib179+LLL79EXV0d+vfvj7i4OCxatAjArWP69NNP4+DBg9ixYweMjIzg5uaG1NRUrX8mqPNQSK25AomIiO6Zj48P+vfvj08++UTuUkgwDHciIhmoVCrY2toiLy9PZ0qF6F4x3ImIiATDq+WJiIgEw3AnIiISDMOdiIhIMAx3IiIiwXT4fe5qtRqXLl2CUqmU5SMxiYiI9JUkSaiqqkLv3r3v+CFYHR7uly5d0vnCCiIiImq+kpKSO34BUYeHe8P3apeUlOh8TSIRERE1TaVSwcHBQZOlTenwcG8YirewsGC4ExERtcLdprV5QR0REZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCMZK7AAIKCwtRVVUldxmdjlKphKurq9xlEBG1OYa7zAoLCzFw4EC5y2gV+24KzPYyxqbvr+NytSR3Oa1SUFDAgCci4TDcZdZwxr5z5064u7vLXE3LmFUUwD1nNqa+sQ3XuuvXPyinT59GeHg4R0yISEgM9/uEu7s7HnzwQbnLaJlLBkAO4O7mBvQeJnc1RET0/3hBHRERkWAY7kRERIIRItxra2tx8uRJ1NbWyl0KEQmK7zOkT4QI9zNnzsDLywtnzpyRuxQiEpS+vs8sX74cCoVC87N8+XK5S2qW8vJyODk5oVu3bnByckJ5ebncJTVbfX09Dh8+jA8//BCHDx9GfX19h9fQ4nDPycnBhAkT0Lt3bygUCmRkZLRDWUREdK8UCgViY2O12mJjY6FQKOQpqJm6d+8OOzs7FBcXo6amBsXFxbCzs0P37t3lLu2u0tPT4eLigkceeQTTp0/HI488AhcXF6Snp3doHS0O95qaGgwdOhTr169vj3qIiKgN3B7glpaWd1x+v+jevTsqKysBAB4eHsjMzISHhwcAoLKy8r4O+PT0dEyePBmenp44evQoqqqqcPToUXh6emLy5MkdGvAtDvegoCDExcXhqaeeao96iIjoHv116H3r1q2QJAkVFRWQJAlbt25tdL37QXl5uSbYKysrkZ+fj/HjxyM/P1+r/X4coq+vr8eCBQvw5JNPIiMjAyNHjkS3bt0wcuRIZGRk4Mknn8TChQs7bIi+3e9zr6urQ11dneaxSqVq8z6uXbsG4NYHk+ibhpob9oE6hj6/Zkge+vS3+teh+BdeeEFr2QsvvIAXX3xRs15MTExHlnZHPj4+AG6dsVtYWGgts7CwgLu7O06fPg0fHx8UFRXJUWKTcnNzUVxcjA8//BAGBtrnzQYGBnj99dfh6+uL3NxcBAQEtHs97R7u8fHx7f7fYXFxMQAgPDy8XftpT8XFxRg1apTcZXQaIrxmSB769Ld6+1B8g27duqG6urqDq7m7hjPyxMTERpevWLECkyZNui/P3H///XcAwODBgxtd3tDesF57a/dwf/311zF//nzNY5VKBQcHhzbtw9HREYB+foRrw8egNuwDdQx9fs2QPPTxb7VhKPt292OwA4CtrS1qamrw2muvYfz48TrLlyxZolnvftOrVy8AQH5+PkaOHKmzPD8/X2u99tbu4W5iYgITE5N27cPMzAyAnn6E6/9r2AfqGCK8Zkge+vC3Ghsbqxmaf//997WG5t9//32t9e4nx44dg52dHU6dOgWVSqU1NK9SqTRTI8eOHZOrxCY9/PDDcHR0xMqVK5GRkaE1NK9WqxEfHw8nJyc8/PDDHVKPEPe5ExHRn/46j/7iiy9CoVBAqVRCoVBo5ttvX+9+YGtrq5lKsLS0xKBBg7Bnzx4MGjRIq/1+PHM3NDTE6tWrkZmZiZCQEK2r5UNCQpCZmYm3334bhoaGHVJPi8/cq6urcfbsWc3joqIi5OXlwcrKCv369WvT4oiIqHUkSdK63e32oXhJuj+/prmiokJzO9zp06cxadIkzTJLS0tUVFTIV9xdTJo0CWlpaViwYAF8fX017U5OTkhLS9Pal/bW4jP3EydOYPjw4Rg+fDgAYP78+Rg+fDjeeOONNi+OiIhaT5KkRj/E5n4N9gYVFRUoKyuDo6MjunbtCkdHR5SVld3Xwd5g0qRJOHv2LLKzs7F7925kZ2ejsLCwQ4MdaMWZe0BAwH3/wiAioltiYmLuu+H35rC1tb3vbndrLkNDww653e1OhJhzd3Nzw/fffw83Nze5SyEiQfF9hvRJu18t3xHMzc15xTMRtSu+z5A+EeLMnYiIiP7EcCciIhKMEMPy+qy2thYAcPLkSZkraTmzigK4Azh95gyuXVbLXU6L8DPliUhkDHeZnTlzBgDw0ksvyVxJy9l3U2C2lzE2rZ6Oy9X6eQeFUqmUuwQiojbHcJdZSEgIgFtX4pqbm8tbTCsFy11AKymVSri6uspdBhFRm1NIHXzTukqlgqWlJSorK3W+0o+IiIia1twM5QV1REREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgmG4ExERCYbhTkREJBiGOxERkWAY7kRERIJhuBMREQmG4U5ERCQYhjsREZFgGO5ERESCYbgTEREJhuFOREQkGIY7ERGRYBjuREREgjHq6A4lSQIAqFSqju6aiIhIrzVkZ0OWNqXDw72qqgoA4ODg0NFdExERCaGqqgqWlpZNLldId4v/NqZWq3Hp0iUolUooFIo22aZKpYKDgwNKSkpgYWHRJtuku+NxlwePuzx43Dsej7kuSZJQVVWF3r17w8Cg6Zn1Dj9zNzAwQN++fdtl2xYWFnwByIDHXR487vLgce94POba7nTG3oAX1BEREQmG4U5ERCQYIcLdxMQEMTExMDExkbuUToXHXR487vLgce94POat1+EX1BEREVH7EuLMnYiIiP7EcCciIhIMw52IiEgwDHciIiLB6H24r1+/Ho6OjjA1NcVDDz2EY8eOyV2S0OLj4zFixAgolUrY2dkhJCQEv/76q9xldToJCQlQKBSIioqSuxThXbx4EeHh4bC2toaZmRk8PT1x4sQJucsSWn19PZYtWwYnJyeYmZnB2dkZb7755l0/T53+pNfh/vHHH2P+/PmIiYnByZMnMXToUIwbNw5lZWVylyasr7/+GnPmzMF3332HAwcO4MaNG3jsscdQU1Mjd2mdxvHjx7Fp0yYMGTJE7lKEd/XqVYwaNQpdunTB/v378csvv2D16tXo0aOH3KUJLTExEe+99x7WrVuH06dPIzExEW+99RaSkpLkLk1v6PWtcA899BBGjBiBdevWAbj1ufUODg545ZVXEB0dLXN1nUN5eTns7Ozw9ddfw8/PT+5yhFddXY0HH3wQGzZsQFxcHIYNG4a1a9fKXZawoqOj8c033yA3N1fuUjqVJ598Ej179sTWrVs1baGhoTAzM8POnTtlrEx/6O2Z+/Xr1/H9999j7NixmjYDAwOMHTsWR48elbGyzqWyshIAYGVlJXMlncOcOXMwfvx4rdc9tZ+9e/fC29sbYWFhsLOzw/Dhw7F582a5yxKer68vDh48iIKCAgDAjz/+iCNHjiAoKEjmyvRHh39xTFu5cuUK6uvr0bNnT632nj174syZMzJV1bmo1WpERUVh1KhRGDx4sNzlCO+jjz7CyZMncfz4cblL6TTOnz+P9957D/Pnz8fixYtx/PhxREZGwtjYGDNmzJC7PGFFR0dDpVLBzc0NhoaGqK+vx4oVK/DMM8/IXZre0NtwJ/nNmTMH+fn5OHLkiNylCK+kpATz5s3DgQMHYGpqKnc5nYZarYa3tzdWrlwJABg+fDjy8/OxceNGhns7Sk1Nxa5du7B79254eHggLy8PUVFR6N27N497M+ltuNvY2MDQ0BClpaVa7aWlpbC3t5epqs5j7ty5yMzMRE5OTrt9hS/96fvvv0dZWRkefPBBTVt9fT1ycnKwbt061NXVwdDQUMYKxdSrVy8MGjRIq83d3R2ffvqpTBV1DosWLUJ0dDSmTZsGAPD09MSFCxcQHx/PcG8mvZ1zNzY2hpeXFw4ePKhpU6vVOHjwIP72t7/JWJnYJEnC3LlzsWfPHhw6dAhOTk5yl9QpjBkzBj///DPy8vI0P97e3njmmWeQl5fHYG8no0aN0rnVs6CgAP3795epos6htrYWBgba8WRoaAi1Wi1TRfpHb8/cAWD+/PmYMWMGvL294ePjg7Vr16KmpgYRERFylyasOXPmYPfu3fjss8+gVCpx+fJlAIClpSXMzMxkrk5cSqVS57qGrl27wtramtc7tKNXX30Vvr6+WLlyJaZMmYJjx44hOTkZycnJcpcmtAkTJmDFihXo168fPDw88MMPP+Cdd97BCy+8IHdp+kPSc0lJSVK/fv0kY2NjycfHR/ruu+/kLkloABr9SUlJkbu0Tsff31+aN2+e3GUIb9++fdLgwYMlExMTyc3NTUpOTpa7JOGpVCpp3rx5Ur9+/SRTU1NpwIAB0pIlS6S6ujq5S9Mben2fOxEREenS2zl3IiIiahzDnYiISDAMdyIiIsEw3ImIiATDcCciIhIMw52IiEgwDHciIiLBMNyJiIgEw3AnIiISDMOdiIhIMAx3IiIiwTDciYiIBPN/nkTON4bYZGgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(6,2))\n",
    "plt.title(\"Number of verified ranks, \" + method)\n",
    "plt.boxplot(N_verified_all.reshape(-1), vert=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrospective rank verification for top K set\n",
    "### Looks fine!\n",
    "- Had to fix a bit, with wrong def of FWER. \n",
    "- One instance where max FWER is higher than alpha, but it's really not unreasonable to see at least this many errors in 50 samples\n",
    "<!-- - works across K, alpha\n",
    "- possible exception: K=5, alpha=0.2, get FWER of 57% when skip_thresh < 0.3. Hard at the boundary here; only rejects 14/50 times. Empirical FWER is very high, 0.38. So it may just be too variable to gauge.\n",
    "    - Kind of a bummer for analysis: \"Hey our method works unless it's a boundary case, in which case it should be taken with a grain of salt... but there's no way of knowing whether it's a boundary case without rerunning a bunch of times\" -->\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alpha: 0.05 \tAverage FWER (out of 30): 0.005 \tMax FWER: 0.080\n"
     ]
    }
   ],
   "source": [
    "importlib.reload(top_k)\n",
    "importlib.reload(helper)\n",
    "K = 5\n",
    "alphas = [0.05, 0.1, 0.2]\n",
    "# skip_thresh = 0.2\n",
    "for alpha in alphas:\n",
    "    fwers_all = []\n",
    "    for i in range(N_pts):\n",
    "        num_false_rejections = 0\n",
    "        true_top_K_set = np.sort(avg_ranks[i,:K])\n",
    "        for j in range(N_runs):\n",
    "            ss_vals, ss_vars = shap_vals[i,j,:], shap_vars[i,j,:]\n",
    "            result = helper.test_top_k_set(ss_vals, ss_vars, K=K, alpha=alpha, abs=True)\n",
    "            if result==\"reject\":\n",
    "                est_top_K_set = np.sort(all_ranks[i,j,:K])\n",
    "                if not np.array_equal(true_top_K_set, est_top_K_set):\n",
    "                    num_false_rejections += 1\n",
    "        fwer = num_false_rejections/N_runs\n",
    "        fwers_all.append(fwer)\n",
    "    print(f\"Alpha: {alpha} \\tAverage FWER (out of {len(fwers_all)}): {np.mean(fwers_all):.3f} \\tMax FWER: {np.max(fwers_all):.3f}\")\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.240\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import binom\n",
    "prob_get_that_max = 1-binom.cdf((50*0.08-1), 50, 0.05)\n",
    "print(f\"{prob_get_that_max:.3f}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# helper.calc_fwer(np.sort(all_ranks[i,:,:K], axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Top-K\n",
    "## Rank stability\n",
    "### from rankshap_vs_ss.py\n",
    "<!-- - alpha = 0.2: Always controls FWER for K=2 and K=5, didn't let finish, fairly conservative\n",
    "- alpha = 0.1: Generally good. Didn't let K=5 finish. 1x on each K, get  FWER > alpha, but by a margin that's not unlikely.  -->\n",
    "\n",
    "- Just changed so it saves the iterations that fail to converge, and computes FWER based on both.\n",
    "- haven't run for alpha=0.2 or 0.1 w/ K=5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['rankshap', 'rankshap_rejection_idx', 'ss_adaptive', 'rankshap_n_samples', 'x_indices'])\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.1\n",
    "top_k_path = join(results_path, \"Top_K\", \"rank\", \"alpha_\" + str(alpha))\n",
    "K = 2\n",
    "fname = \"rankshap_vs_ss_k\" + str(K)\n",
    "\n",
    "with open(join(top_k_path, fname), 'rb') as f:\n",
    "    results = pickle.load(f)\n",
    "print(results.keys())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6\n",
      "[0.0, 0.0, 0.0, 0.0, 0.04, 0.0]\n"
     ]
    }
   ],
   "source": [
    "top_K_rs = results['rankshap']\n",
    "rejection_idx = results['rankshap_rejection_idx'] # Doesn't exist yet\n",
    "\n",
    "print(len(top_K_rs))\n",
    "print([helper.calc_fwer(top_K, digits=3, rejection_idx=rejection_idx[i]).item() for i, top_K in enumerate(top_K_rs)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### top_k_experiment.py\n",
    "<!-- OUTDATED:\n",
    "- RankSHAP alpha=0.1: Always controls FWER, not too conservative, made it to 30, and fixed saving error with FWERs\n",
    "- RankSHAP alpha=0.2: Good but didn't make it to 30\n",
    "- SPRT-SHAP alpha=0.1: Good but didn't make it to 30; conservative -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.  , 0.  , 0.  , 0.  , 0.02, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,\n",
       "       0.  ])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = 0.1\n",
    "top_k_path = join(results_path, \"Top_K\", \"rank\", \"alpha_\" + str(alpha))\n",
    "K = 2\n",
    "# fname = \"rankshap_census_K\" + str(K)\n",
    "fname = \"sprtshap_census_K\" + str(K)\n",
    "\n",
    "with open(join(top_k_path, fname), 'rb') as f:\n",
    "    results = pickle.load(f)\n",
    "\n",
    "# top_K_all = results['ranks']\n",
    "# idx = results['rejection_idx'] #\n",
    "# fwers = np.array([helper.calc_fwer(top_K, digits=3, rejection_idx=idx[i]).item() for i, top_K in enumerate(top_K_all)])\n",
    "results['fwers']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('hello')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fwers = results['fwers']\n",
    "\n",
    "# fwers = results['fwers'].item()\n",
    "# fwers = np.array(list(fwers.values()))\n",
    "# results['fwers'] = fwers\n",
    "# with open(join(top_k_path, fname), 'wb') as f:\n",
    "#     pickle.dump(results, f)\n",
    "# fwers\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rankings",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}