notebooks/06_fig6_cd_seg_plots.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import os\n",
    "from os.path import join as oj\n",
    "import sys\n",
    "sys.path.append('../src')\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import torch\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn import metrics\n",
    "from sklearn.linear_model import RidgeCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "# plt.style.use('dark_background')\n",
    "import data\n",
    "from skorch.callbacks import Checkpoint\n",
    "from skorch import NeuralNetRegressor\n",
    "from config import *\n",
    "from tqdm import tqdm\n",
    "import pickle as pkl\n",
    "import train_reg\n",
    "from math import floor\n",
    "from copy import deepcopy\n",
    "import config\n",
    "import models\n",
    "import pandas as pd\n",
    "import features\n",
    "from scipy.stats import skew, pearsonr\n",
    "import outcomes\n",
    "import neural_networks\n",
    "from sklearn.model_selection import KFold\n",
    "from torch import nn, optim\n",
    "from torch.nn import functional as F\n",
    "import interpret\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "import viz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# prepare data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1/1 [00:15<00:00, 15.72s/it]\n"
     ]
    }
   ],
   "source": [
    "# currently trained only on 'clath_aux+gak_a7d2_new'\n",
    "dsets = ['clath_aux_dynamin']\n",
    "splits = ['test']\n",
    "feat_names = ['X_same_length_normalized'] # + data.select_final_feats(data.get_feature_names(df))\n",
    "              #['mean_total_displacement', 'mean_square_displacement', 'lifetime']\n",
    "meta = ['cell_num', 'Y_sig_mean', 'Y_sig_mean_normalized']\n",
    "dfs, feat_names = data.load_dfs_for_lstm(dsets=dsets, splits=splits, meta=meta, normalize=False)\n",
    "\n",
    "# load model\n",
    "p = 1\n",
    "results = pkl.load(open(config.FINAL_MODEL, 'rb'))\n",
    "dnn = neural_networks.neural_net_sklearn(D_in=40, H=20, p=p-1, arch='lstm')\n",
    "dnn.model.load_state_dict(results['model_state_dict'])\n",
    "\n",
    "# load data\n",
    "# df = dfs[('clath_aux+gak_a7d2', 'train')]\n",
    "df = dfs[('clath_aux_dynamin', 'test')]\n",
    "X = df[feat_names[:p]]\n",
    "y = df['y_consec_thresh']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# segmentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 55/8451 [00:00<00:15, 549.09it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "acc 0.8409655662051828\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 8451/8451 [00:15<00:00, 536.49it/s]\n",
      "<ipython-input-3-0949b7040d79>:48: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  all_res = np.array(all_res)\n"
     ]
    }
   ],
   "source": [
    "track_name = 'X_same_length_normalized'\n",
    "track_name_unnormalized = 'X_same_length'\n",
    "df = dfs[('clath_aux_dynamin', 'test')]\n",
    "df['preds'] = dnn.predict(df[feat_names[:1]])\n",
    "df = df.sort_values(by='preds')\n",
    "x = df[feat_names[:1]]\n",
    "y = df['y_consec_thresh'].values\n",
    "preds = dnn.predict(x)\n",
    "n = df.shape[0]\n",
    "\n",
    "m = 40\n",
    "xtrack_unnormalized = df[track_name_unnormalized]\n",
    "xtrack = x[track_name]\n",
    "xtrack_t = torch.tensor(np.array(list(xtrack.values)), dtype=torch.float)\n",
    "xfeats = x[[c for c in x.columns if c != track_name]]\n",
    "xfeats_t = torch.tensor(np.array(xfeats).astype(float), dtype=torch.float)\n",
    "\n",
    "\n",
    "# sort things\n",
    "print('acc', np.mean((preds > 0) == y))\n",
    "\n",
    "# calculate all cd scores\n",
    "cd_score_path = f'../data/outputs/all_cd_scores_clath_aux_dynamin.pkl'\n",
    "if os.path.exists(cd_score_path):\n",
    "    all_cd_scores = pd.read_pickle(cd_score_path)\n",
    "else:\n",
    "    all_cd_scores = {}\n",
    "    for s in tqdm(range(m)):\n",
    "        for e in range(s+1, m+1):\n",
    "            all_cd_scores[(s, e)] = interpret.calc_cd_score(xtrack_t, xfeats_t, s, e, dnn.model)\n",
    "    with open(cd_score_path, 'wb') as handle:\n",
    "        pkl.dump(all_cd_scores, handle, protocol=pkl.HIGHEST_PROTOCOL)\n",
    "        \n",
    "# find best segmentations\n",
    "all_res = []\n",
    "for i in tqdm(range(len(x))):\n",
    "    cd_scores_list = np.zeros((m, m))\n",
    "    for s in (range(m)):\n",
    "        for e in range(s + 1, m + 1):            \n",
    "            cd_scores_list[s][e-1] = all_cd_scores[(s, e)][i]\n",
    "    res, paths = interpret.max_abs_sum_seg(cd_scores_list, min_length=5)\n",
    "    lt = np.sum(np.array(df[track_name_unnormalized].values[i]) != 0)\n",
    "    #all_res.append(paths[lt-1])\n",
    "    all_res.append(paths[m - 1])\n",
    "    #print(res[lt-1], paths[lt-1])\n",
    "\n",
    "# paths is a list of lists, where each list contains all the starting point for segments    \n",
    "all_res = np.array(all_res)    \n",
    "    \n",
    "def extract_segs_and_scores(all_res, all_cd_scores, i):\n",
    "    '''\n",
    "    Params\n",
    "    ------\n",
    "    all_res: array_like\n",
    "        list of list of starting points for all segments\n",
    "    '''\n",
    "    segs = [(all_res[i][j], all_res[i][j+1]) # converting starting points to list of (start, end) pairs\n",
    "            for j in range(len(all_res[i]) - 1)\n",
    "           ] + [(all_res[i][-1], 40)]\n",
    "    scores = [all_cd_scores[(s, e)][i][0]\n",
    "              for (s, e) in segs]\n",
    "    return segs, scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for which_tracks in ['aux-', 'aux+', 'random']: # random, neg, pos\n",
    "#     C, R = 3, 5\n",
    "    for setup in ['small', 'big']:\n",
    "        if setup == 'small':\n",
    "            R = 3\n",
    "            C = 2\n",
    "            plt.figure(dpi=150, figsize=(2.6 * C, 1 * R), facecolor='w')\n",
    "        else:\n",
    "            R = 5\n",
    "            C = 3\n",
    "            plt.figure(dpi=150, figsize=(1.8 * C, 1 * R), facecolor='w')\n",
    "        # calculate vabs\n",
    "        scores_all = [extract_segs_and_scores(all_res, all_cd_scores, i)[1] for i in range(all_res.size)]\n",
    "        vmax = max([max(x) for x in scores_all])\n",
    "        vmin = min([min(x) for x in scores_all])\n",
    "        vabs = max(abs(vmax), abs(vmin)) / 3 # note this thresholds things\n",
    "\n",
    "        \n",
    "        i = 0\n",
    "        for r in range(0, R):\n",
    "            for c in range(0, C):\n",
    "                ax = plt.subplot(R, C, i + 1)\n",
    "                if which_tracks == 'random':\n",
    "                    arg = np.arange(n)[n//2:] #np.random.permutation(np.arange(n))\n",
    "                elif which_tracks == 'aux-':\n",
    "                    arg = np.arange(n)\n",
    "                elif which_tracks == 'aux+':\n",
    "                    arg = np.arange(n)[::-1]\n",
    "                segs, scores = extract_segs_and_scores(all_res, all_cd_scores, arg[i])\n",
    "            #     print('scores', segs, scores, preds[arg][i])\n",
    "                interpret.plot_segs(track_segs=segs,\n",
    "                                    cd_scores=scores,\n",
    "                                    xtrack=xtrack_unnormalized.iloc[arg[i]],\n",
    "            #                         pred=preds[arg][i],\n",
    "            #                         y=y[arg][i],\n",
    "                                    vabs=vabs,\n",
    "                                    cbar=False,\n",
    "                                    xticks=(c==0 and r==R-1),\n",
    "                                    yticks=i==0)\n",
    "\n",
    "    #             plt.subplot(R, C, 1)\n",
    "    #             plt.ylabel('Fitted amplitude')\n",
    "                plt.subplot(R, C, 1 + C * (R - 1))\n",
    "                plt.xlabel('Time')\n",
    "\n",
    "                if which_tracks in ['aux+']:\n",
    "                    plt.ylim((-100, 25000))\n",
    "                elif which_tracks in ['random']:\n",
    "                    plt.ylim((-100, 10000))\n",
    "                else:\n",
    "                    plt.ylim((-100, 1000))\n",
    "\n",
    "                '''\n",
    "                plt.subplot(R, C, 2)\n",
    "                plt.title(which_tracks\n",
    "                             .replace('random', 'Random')\n",
    "                             .replace('aux-', 'Abortive')\n",
    "                             .replace('aux+', 'Valid') + \\\n",
    "                             ' events')\n",
    "                '''\n",
    "                plt.suptitle(which_tracks\n",
    "                             .replace('random', 'Predicted as unsure')\n",
    "                             .replace('aux-', 'Predicted as most abortive')\n",
    "                             .replace('aux+', 'Predicted as most succesful') + \\\n",
    "                             ' events', y=.997)\n",
    "                plt.figtext(s={\n",
    "                    'random': 'C',\n",
    "                    'aux-': 'A',\n",
    "                    'aux+': 'B'\n",
    "                }[which_tracks], x=0.05, y=.997, fontweight='bold', fontsize=12)\n",
    "                i += 1\n",
    "#         plt.tight_layout()\n",
    "#         plt.tight_layout(rect=[0, 0.03, 1, 0.5]) \n",
    "#         plt.subplots_adjust(top=0.85)\n",
    "    \n",
    "\n",
    "        viz.savefig('interpretations_' + which_tracks + f'_{setup}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "../src/viz.py:51: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n",
      "  plt.savefig(oj(DIR_FIGS, 'fig_' + s + '.pdf'), bbox_inches='tight')\n",
      "../src/viz.py:51: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n",
      "  plt.savefig(oj(DIR_FIGS, 'fig_' + s + '.pdf'), bbox_inches='tight')\n",
      "/system/linux/anaconda3.8/lib/python3.8/site-packages/IPython/core/pylabtools.py:132: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n",
      "  fig.canvas.print_figure(bytes_io, **kw)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAABUCAYAAAChrxGdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAUZklEQVR4nO2debxdVXXHv78AAiESMICRhpoPAgawgxVQizIIMiiICSFoERtlEqfWailakaC2HxUDlioaQKYKhMGgMs/BEZFB0aJBJalGQiSDTBkA3+ofa5+88847995z7zv3vpe89f18zufes88e19l7r332KDMjCIIgCOpizHBHIAiCINiwCMUSBEEQ1EooliAIgqBWQrEEQRAEtRKKJQiCIKiVUCxBEARBrWzcjuUHZuxhSKCkjyS/1v1vbq6iOfizAXa7Z65O/CmmjfbSXPRb3UhrG3FRx3Efovmg8Kumud64DMqDJXFsJX9DZNP0zSCbsN9nRjZ732Cdnb5kj2TWB/1uM3+Avpwdy9npy9npDzP5bTm/6Q+/L2fH/Rhsxwz6GJwOM8vZZ0Ba+3J28mktS38+3nWlP5/Wvlw6WqW1u+kvj0tfwc91cSn4mZdf5bg0kFdfQab96R8o09K82ST/9Me93wzg/k8dnQrIQOKLJQiCIKiVUCxBEARBrYRiCYIgCGolFEsQBEFQK6FYgiAIgloJxRIEQRDUSiiWIAiCoFZCsQRBEAS1EoolCIIgqJVQLEEQBEGthGIJgiAIaiUUSxAEQVAroViCIAiCWgnFEgRBENRKKJYgCIKgVkKxBEEQBLUSiiUIgiColVAsQRAEQa2EYgmCIAhqJRRLEARBUCuhWIIgCIJaCcUSBEEQ1EooliAIgqBWQrEEQRAEtRKKJQiCIKiVUCxBEARBrYRiCYIgCGolFEsQBEFQKzKz7gcinWhm53U9oFFOyLk3hJx7Q8i5N3RDzr36YjmxR+GMdkLOvSHk3BtCzr2hdjlHV1gQBEFQK6FYgiAIglrplWKJftLeEHLuDSHn3hBy7g21y7kng/dBEATB6CG6woIgCIJa6YpikXS0pHmSlkgySTPbcLu3pB9LWiNpoaQPdyOOGwqSTpD06ySv+yUdUMHNrPReitchvYjzSEXSbpLukLRK0mOSPi1powruxku6SNJKSU9KukzShF7EeX2kEzlLmtwgz87tVbzXNyTtJGmOpIck/VnS/IruhpyfN+4oxq2ZDkwGrgeOr+pI0k7ALcndx4G9gLMkrTKzC7oQz/UaSe8EvgbMAr4PvAe4XtKeZvaLFs6fBIqK5Je1R3I9QdLWwO3Aw8ARwCuA2Xjj65MtnF8F7ILn9T7g88C3gDd2K77rK0OUM8DHgB/k7pfVHccNiN2BtwD3AJu04W7o+dnMar+AMel3HGDAzIru5gCPABvnzM4Ffk8aD4prgLwWABfm5Q78HPhGC3ezgGXDHf+RdOENmZXAljmzU4BVebMSd69PeXyfnNleyezA4U7XSLuGIOfJSaaHDXca1pcrq4fT/2uA+RXc1JKfu9IVZmZ9HTo9FJhnZi/kzOYCk4BXDTliGxCSdsRbFVdlZknuV+NyDNrjUOAWM3sqZzYX2BzYt4W7pWb23czAzO4FFhLvoYxO5Ry0SYf1cC35ecQM3kvaAtgB+FXhUdY9M6W3MRrxZPIok9dLJG3bwv1WkpZJel7Sg5Km1R/F9YopFGRpZr/DW9LN8t4gd4lftnA3WulUzhkXpfGCJZLOkrR5NyI5iqklP48YxQJslX7/VDBfmX637mFc1gcyeXQir9/g3Q9HAUcCjwHfHOXKZWsGyxJcns1k2am70Uqn8loLfAU4DjgA7zY/Gf/aCeqjlvxcafBe0njgZa3smVmZpgsq0is5m9k3CuFeB/wQ+BQwbyh+B0E3MLMlwAdzRvMlLQXOlfQ3ZvazYYpaUELVWWFHAedXsKchxCXTkuML5pmWXMmGTztyzuQxnoEtjLblZWYmaR7weUkbmdmfq7rdgFjJ4LwHLs9mslwJlHU7tnI3WulUzmVcg0/ueQ0QiqUeasnPlbrCzOwCM1Orq2qgDcJ4Fp/9VezHazSWsMHRppwzeZTJa4WZPdFu8OkarfyKgiwl7QCMpXneG+Qu0aiverTTqZzLsMJvMHRqyc8jaYwF4CZgamGx1NG4wmm1LmNUYWaP4lOzj8rMJI1J9ze145ck4WMtPxulXyvgMjtY0otzZkcDq4G7W7ibKOkNmYGkPYAdafM9jBI6lXMZ09Pv/XVELADqys9dmj+9G/7S34W3Jr6c7vfN2dkXeKFgthPwDHA5sD8+wPw8cPxwzwkfiRfwTuDP+MKy/YGL8QL6qhZyvhv4MHAQMBW4EV8I9bbhTtMwynJrYAlwG3AgfkbFM8BnC/Z+A3y9YHYL8CgwDXg7vr7oe8OdppF4dSpnfO3V7CTjA4FPp7z+zeFO00i98K/A6en6EfC/ufuxZXJOZkPOz91K0Cz6u1by1/ycnf2S2X4Ft28A7gXWAIuADw/3CxrJF3BCyhxrgQeAAwrPB8kZ+HrKOKuBZ4HvAYcOd1qG+8IbRHcmuSwBPgNsVLCzCLi4YLYVcBE+1vUU3jDaZrjTM1KvTuQMvAO4D98x4rmU5z8NbDrc6RmpF/2LSsuuyWVyTmZDzs+xu3EQBEFQKyNtjCUIgiBYzwnFEgRBENRKKJYgCIKgVkKxBEEQBLUSiiUIgiColVAsQRAEQa30XLFIOlLSnZL+JGmtpEfS9tfbdyGsXdIxvFu1tj3g+NPDcmaLJH2xpvicKOntJea1hVEhDi9KMvnbXoTXTdK2/7Ny9/MlXdOG+9L8IWlmygfjaoxukEjHjz+QjtOuvN6hnfci6RxJF5WYt6x/Ckcfr5b0O/lR64dXT+XIRdJ0SQuqHLvdKT1VLJJm4wdTPQoci6/8PhvfBvsrXQhyF+B0+rfkb8US/AS173chLuCrjAcpFnz1+zldCrPIi3CZrPeKpYT34ycUVqVR/rgBzweraopXMJA5+OK7g3E510rae+wE/EjdvHk79c/sFLeDgFPxRZnflnRh3fEdBubhG9ke260AunXm/SCStv8X4Dgzy7+cuyWdh7/AYUPSZma2Bj8fuqeY2YO9DnO4SfuTbZpkXgtm9nBN/jwBtLuJZ1CdKcB5Ztbu3mBVeR/wgOWOl+ig/llkZvm64HJJtwJfl3S3mV3Spbh3RDvlycz6JF0KfAjfBqp+eri9wJ3A/RXtbgNcAizHW43zgT0KdhYBXwQ+AizGt3SeC2yVnu/H4G0MFqVnM9P9Xsnv1cBplJyrnQvnNOBxfF+jy4DxOTuZf+PK4pj+zy+Jz8yivZzbGfj59WvxTTj/A9i4JMy/wvddehbffXRaC9k2296hpdybbBvxD8D/AE8DfwROL9ibBSzDt+z5SUrXsenZG/H9y1alsM8HXlxwvw++NfoafNPBv0/+zcrZmQ9cU3D318B1eAv5GXy7oDdXzB/jCmmcgbe2n8Tz3BnkzhVPdo8Cfo3nqbuAV+ffdRM5fhzfpmQNsBS4GZhYNX/lzKamNK5OsrwReHkreeSevwQ4L8VhDX5Oz2sLYRwHPJzCWJbe3e6t0tJA5hfn8uUHy/JMq3JWIsvfAx8ZQv0zKC65Z/cA97Rw/zY8jz6L10s/ZuBefRslGT2Cl4PFDN5W5YMpH61NsiymZxZDK0+7pnTu3iwtnV496QqTtAleEdxc0cm38M/kj+E7n44B7pK0U8HeDPwz9kTg34DDgP9Mzx5I7sE3U3s9XujyXIEXsrcA1zeJzzvxje9OwFs9bwUuqJiWjPfjFf+NKS6vx7tcBiHpIODKlIYjgP9OaflyifXLge/gafs1MFfSpCbxeFP6/WwuHkuSWVW5l3EmnpGn4xn5dEkfKNgZiyuuC4BDgHsl7Q3cjivt6cA/4+9jXf946v++CViR7MzBlfvYZhGSNAX4AX542vtwGV2LH4FdJX8U+QJeGU8HvoEfjJbtsJvtAjs3+T0Vfy9XtvATSe8GPgGchcv/ZLwy2aKV24I/x+LdHL/Fy8Z78Mpr2/S8mTyQtCn+Lg4E/hXvtn0CuF3SxGRnH+BreCPiUOC9uPIZXyEtD9Df9ZV1NX2mnTRWkMErgUkpTplZu/VPM24DXpP8LAv/Ffg5MXcChwPH4HXLS3LW5uCNkqvwOuuj5PKypBPwMv+d5MfVwGxJpxaC66g8AZjZL3Gld2Bbqa9KN7RViQafiGvHkyrYPSTZzWv4LfAMPqfQWvstA1vxXwIez90fRq5FXtLy+aeC+WTKv1hWkGsl4ZmlD9i1WUuKQosS30Tv4pI0F+3dA9xVsHMKvpPxpEKY783ZmYDvZPy+JvIdR0kLuqrcS/zLZHZrwfx84A+kFj39G5MeUbD3vZK0vinZfVW6/wLe8hpbeAdGky8WvOGwGNi8Qdxb5Y/iF8ulBXs/Bebm7q/Gj3dQ4b01/WLBGwwNd+mtkr/wRsAfgHlN/Gklj+PwsYSdc2Yb4+XszHT/MZq0/FulJdkp+zqp5YsF/3I2YIucWeX6p1Fccs9OSs9f2uD5dGB5E7+nJPelm+vm3uNFBfNz8S/lzYZangrl5bIqMmn36vWsMKtgZy/gj5brfzU/BOx6/LMvz11m9kLu/mFgu0atiRJKvxhKuM3MnsndX4sPfu1Z0X1l0kyNv8MrqTxX4pmuONh5a/bHzJbj3VDNvlga0Y7cy7i2cD8P2L4QFyN3poOksXh6rpK0cXbhkyeex08GzOJ2m5nlB9OL4ZXxJuBKM1tdwW4Vbi3cP8zA9O0JXGep1Ca+U8HfnwJvkXSGpL06nK3zSlzeFzWx00oeB+JdOAtz7wK8W2WPXFxfLelsSftIelEX0jIUJgJrUt4tUqX+aUWrAw1/DoyXdImkgyQVvzr3T78XN3A/CX+PZeV/S7zrO6PT8pSxDJdX7fRKsSzH+wD/soLdl+GVY5GlDPychIFH8oK3tgRsWjFeSyvaGxCfVME9Q4Xz6TtgG2ATBsctu68ig806CLcduZdRdJvd52W00syey91vjfc3n4tn/Oxai8tgh2RvYtH/3DtoxgT6u/nqoJWsJzJ40L/KJIAL8e6jGXh//FJJn22zUp6Qfpult5U8tgFex8B38TzepbYDgJndnu73wVu8yyR9JVeB1pGWobAZnn/ytFP/tOIvcJmsKHtoZgvw7usd8W7vZZIul5Qd9zsBeNbMnmrgf1ZeqpT/TstTxlo6qyta0pNZYWb2vKQf4H2un2xhfQmwXYn5S2nwMocStYr2BsQntQzG0V9Is5kYxdbb1rTPMjwzFGXw0vRbtwwyhir3otvsPl+RFeX9p2Q2Cy+ERR5Lv48X/c+9g2YspzvKvxGPM/i88LLzwwdgZn34tNez01TZY/DJGovx8Ywq+Wt5+m2W3lbyWIF3155c8mxdZW0+I+qSVFlOS3F/Gji1QloasZZ6ys8KYEtJY1Jc2q1/WnEQ3hX4fCMLZnYDcIOk8fh47JfwMZN34O9gC0lbNlAuWXmpUv47LU8ZW9Gl+qSXXWFfAvaQ9I/FB5LGSDok3f4Y787aJ/d8LP6C2l1fkmnzoWrlNxcWZU3FX+B96X5x+t01syDptfinazE+TeNifjTw/eSOHE7MwMd1ftRWzAfTSCZDlXtx4HsaXkgWl9gF1nW13QO80szuK7mygvAT/B3kB+tbDbQD3AHMkNRI5nXlj4yfAIenqZ8Zb2vHAzP7vZl9Dh/w3i0ZV8lfC/C++UHlK0credyBn+L6u5J38fOSuD5hZnPwfv3dSp6XpaURixmYvjH4xJx2WYD3Wry8YF61/mmIpPfi3bJfrRIRM3vSzC7Hu22z9N+Zft/dwNliXAGUlf+n8K62RuFVLU8Zk/HJHbXTs3UsZnadpLPweeB7A9/GuzKm4DNUFgE3m9ktkn4IXJlmQSzHBww3x2cetcOC9HuSpLnAqrICUoHVeAvkTLzFdyZwrfWvm7gXL9TnSDoN/1w9Bc8IeX6Fn/d9MJ6uhWlcpMjpwC3ylcNz8X7VzwDnm1nDiroKZvacpIV4BfMLvDX8UA1y313SHOCbeDfJcfjkiL4W7k4B7pDUh8+meRrvsngr8O9m9gheKXwAuD7loe3x6Zqtxk7OwCv776bFccvx6b/Lzdcy1JU/Mj6PK+i56d3tis8kBG8UlJLktgKvFJ7E++F3xmc6QoX8Zb424RTgMkmX4QP1ho+rXGFm91WQx6V4WZwv3wniUbzrZi98UszZks5I4c/Hv65fjR9/fWrFtDTiWuADkh5M4R7P4IZZFe7FJ7C8BliYk0+l+ifnz2RJr8O7kCbh3VszgAvN7NJGgUs6CR/nuBlXEDvjSuLSFI8Fad3MbEnbAd/Fvxymm9k70nucBcyRtByfhbYv/hX5CWu9TqVKeSJ1XU7Bl1HUTzdmBDS7gCPx+f3ZEaOP4OtEJubsbIu/iJV45XE3sGfBn0UMnsM/k8KsEXwq3//hmW1RI3vJfDLls8Jm45+XS/G56VeQ1svk7O2JF9pVwIPA3sU44v2ut6e0r5sp1CAtR+Otk+fwVkyjdSwt1zaUvIODgIdwpbJuVlQVuZf4lcnsmCSXp/FxhTMYODtqFrkZPgU/XosXxKeSfB/Gp6vm1wrtl+K8Fh8g3pvq61huTPF6Gq/4D8g9b5k/yvJFMr8YuK9gNoP+NRzfxwfEDXh7ExnOxKcBr0j55yF8IV9b+SvZm4Z/8a7BFccNDF7H0kwe44H/wteCZHlvHrB3en4Y/mXzRApjAa5U1EZaymaAjcOnzq7AuxQ/mfJQJ+tYrqdwjnub9Y/lrjVJFvOAwyvUb9kygseS24V4g2PTnJ2N8HGoR3MyvrDgz4dSPnou2StdxzKE8jQ1vf8tWqWpkyuOJg6GhKTJeOE53MyarQUalUh6F77mY0czW9jKfjB0JE3F13Zsb2bFgfwAkHQFPong+G74H7sbB0GNSPqqpGmS9pf0UXzQ9oZQKj3l2/gXQ9f2wlqfSZMqjgA+160wejbGEgSjhAn4dM8JeFfUlXi/d9AjzMcpTsDX9gSDmYQvov5NtwKIrrAgCIKgVqIrLAiCIKiVUCxBEARBrYRiCYIgCGolFEsQBEFQK6FYgiAIgloJxRIEQRDUyv8DPjnI84txun0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_axes([0.05, 0.80, 0.9, 0.1])\n",
    "cb = mpl.colorbar.ColorbarBase(ax, orientation='horizontal', \n",
    "                    norm=matplotlib.colors.Normalize(vmin=-vabs, vmax=vabs),\n",
    "                               cmap=viz.cmap)\n",
    "cb.ax.set_xlabel('Contribution to predicting successful (CD score)')\n",
    "cb.outline.set_visible(False)\n",
    "viz.savefig('interpretations_cbar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**make simple illustrative plot**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 98\n",
    "plt.figure(figsize=(2.4, 1.9), facecolor='w', dpi=300)\n",
    "segs, scores = extract_segs_and_scores(all_res, all_cd_scores, i)\n",
    "cb = interpret.plot_segs(track_segs=segs, cd_scores=scores, xtrack=xtrack_unnormalized.iloc[i])\n",
    "cb.set_ticks([])\n",
    "plt.yticks([])\n",
    "plt.xticks([])\n",
    "plt.show()\n",
    "#               pred=dnn.predict(x)[i],\n",
    "#               y=y.values[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**recalculate with normalized scores (there is some err here)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cd_scores_normalized = {}\n",
    "for s in tqdm(range(m)):\n",
    "    for e in range(s + 1, m + 1):\n",
    "        x = all_cd_scores[(s, e)]\n",
    "        x = x.reshape(1, -1)[0]\n",
    "        all_cd_scores_normalized[(s, e)] = np.abs((x - np.mean(x)) / np.std(x))\n",
    "        \n",
    "for i in range(len(df)):\n",
    "    cd_scores_list = np.zeros((m, m))\n",
    "    s0, e0 = 0, 0\n",
    "    m = 0\n",
    "    lt = np.sum(np.array(df[track_name_unnormalized].values[i]) != 0)\n",
    "    for s in (range(lt)):\n",
    "        for e in range(s+1, lt):            \n",
    "            cd_scores_list[s][e-1] = all_cd_scores_normalized[(s, e)][i]\n",
    "            if cd_scores_list[s][e-1] > m:\n",
    "                m = cd_scores_list[s][e-1]\n",
    "                s0, e0 = s, e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(f'results/all_cd_scores_clath_aux+gak_new_test.pkl', 'wb') as handle:\n",
    "    pkl.dump(all_cd_scores, handle, protocol=pkl.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# interpret one pred at multiple scales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_segs(track_segs, cd_scores, xtrack, pred):\n",
    "    cm = sns.diverging_palette(22, 220, as_cmap=True, center='light')\n",
    "    vabs = np.max(np.abs(cd_scores))\n",
    "    # plt.plot(xtrack, zorder=0, lw=2, color='#111111')\n",
    "    for i in range(len(track_segs)):\n",
    "        (s, e) = track_segs[i]\n",
    "        cd_score = cd_scores[i]\n",
    "        seq_len = e - s\n",
    "        xs = np.arange(s, e)\n",
    "        norm = matplotlib.colors.Normalize(vmin=-vabs, vmax=vabs)\n",
    "        if seq_len > 1:\n",
    "            cd_score = [cd_score] * seq_len\n",
    "            plt.plot(xs, xtrack[s: e], zorder=0, lw=2, color=cm(norm(cd_score[0])), alpha=0.5)\n",
    "        plt.scatter(xs, xtrack[s: e],\n",
    "                    c=cd_score, cmap=cm, vmin=-vabs, vmax=vabs, s=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "track_name = 'X_same_length_normalized'\n",
    "track_name_unnormalized = 'X_same_length'\n",
    "num = 100\n",
    "x = X.iloc[num: num + 1]\n",
    "xtrack_unnormalized = df.iloc[num: num + 1][track_name_unnormalized]\n",
    "xtrack = x[track_name]\n",
    "xtrack_t = torch.tensor(np.array(list(xtrack.values)), dtype=torch.float)\n",
    "xfeats = x[[c for c in x.columns if c != track_name]]\n",
    "xfeats_t = torch.tensor(np.array(xfeats).astype(float), dtype=torch.float)\n",
    "pred = dnn.model(xtrack_t, xfeats_t).item()\n",
    "#print(f'pred {pred:0.2f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "  0%|          | 0/6 [00:00<?, ?it/s]\u001b[A\n",
      " 17%|█▋        | 1/6 [00:02<00:11,  2.35s/it]\u001b[A\n",
      " 33%|███▎      | 2/6 [00:04<00:08,  2.24s/it]\u001b[A\n",
      " 50%|█████     | 3/6 [00:05<00:05,  1.94s/it]\u001b[A\n",
      " 67%|██████▋   | 4/6 [00:06<00:03,  1.58s/it]\u001b[A\n",
      " 83%|████████▎ | 5/6 [00:06<00:01,  1.26s/it]\u001b[A\n",
      "100%|██████████| 6/6 [00:07<00:00,  1.17s/it]\u001b[A\n"
     ]
    }
   ],
   "source": [
    "T = 40 # seq len\n",
    "# track_segs = [(s, s + 1) for s in range(T)]\n",
    "DIV = 40\n",
    "track_segs_fourty = [(floor(s * T / DIV), floor((s+1) * T / DIV))\n",
    "                    for s in range(DIV)]\n",
    "\n",
    "DIV = 20\n",
    "track_segs_twenty = [(floor(s * T / DIV), floor((s+1) * T / DIV))\n",
    "                    for s in range(DIV)]\n",
    "DIV = 10\n",
    "track_segs_tenths = [(floor(s * T / DIV), floor((s+1) * T / DIV))\n",
    "                    for s in range(DIV)]\n",
    "DIV = 6\n",
    "track_segs_fifths = [(floor(s * T / DIV), floor((s+1) * T / DIV))\n",
    "                    for s in range(DIV)]\n",
    "DIV = 4\n",
    "track_segs_quarters = [(floor(s * T / DIV), floor((s+1) * T / DIV))\n",
    "                    for s in range(DIV)]\n",
    "DIV = 1\n",
    "track_segs_full = [(floor(s * T / DIV), floor((s+1) * T / DIV))\n",
    "                    for s in range(DIV)]\n",
    "track_segs_list = [track_segs_fourty, track_segs_twenty, track_segs_tenths,\n",
    "                   track_segs_fifths, track_segs_quarters, track_segs_full]\n",
    "cd_scores_list = [[interpret.calc_cd_score(xtrack_t, xfeats_t, s, e, dnn.model)\n",
    "                  for (s, e) in track_segs]\n",
    "                  for track_segs in tqdm(track_segs_list)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(dpi=200, figsize=(12, 6))\n",
    "R, C = 2, 3\n",
    "for i, track_segs in enumerate(track_segs_list):\n",
    "    ax = plt.subplot(R, C, i + 1)\n",
    "    cd_scores = cd_scores_list[i]\n",
    "    interpret.plot_segs(track_segs, cd_scores, xtrack_unnormalized.iloc[0], pred, y=0)\n",
    "    \n",
    "    if i == C * (R - 1):\n",
    "        plt.xlabel('Time')\n",
    "    else:\n",
    "        plt.xticks([])\n",
    "    if i == 0:\n",
    "        plt.ylabel('Clath Amplitude')\n",
    "    else:\n",
    "        plt.yticks([])\n",
    "    if i  == C - 1:\n",
    "        plt.colorbar(label='CD Score')\n",
    "    else:\n",
    "        plt.colorbar()\n",
    "    if i ==  C * R - 1:\n",
    "        plt.text(0.6, 0.9,\n",
    "                 f'Pred: {pred:.2f}', fontsize='x-large', transform = ax.transAxes)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}