rev_corr_mnist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2628a7ec",
   "metadata": {},
   "source": [
    "# Testing reverse correlation\n",
    "\n",
    "Taken from [https://doi.org/10.1016/j.jneumeth.2021.109297](https://doi.org/10.1016/j.jneumeth.2021.109297)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec9bb8b2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import librairies\n",
    "\n",
    "import struct\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stt\n",
    "\n",
    "from sklearn.model_selection import StratifiedShuffleSplit\n",
    "\n",
    "# classifiers\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# PCA\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "from numba import jit\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "\n",
    "font = {'family' : 'DejaVu Sans',\n",
    "        'weight' : 'regular',\n",
    "        'size'   : 18}\n",
    "mpl.rc('font', **font)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e78fcf9",
   "metadata": {},
   "source": [
    "Load the MNIST dataset: images and labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aef4882d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = './'\n",
    "\n",
    "# use test set as data here\n",
    "with open(data_dir+'t10k-images-idx3-ubyte','rb') as f:\n",
    "    magic, size = struct.unpack(\">II\", f.read(8))\n",
    "    nrows, ncols = struct.unpack(\">II\", f.read(8))\n",
    "    data = np.fromfile(f, dtype=np.dtype(np.uint8).newbyteorder('>'))\n",
    "    data = data.reshape((size, nrows, ncols))\n",
    "\n",
    "with open(data_dir+'t10k-labels-idx1-ubyte','rb') as f:\n",
    "    magic, size = struct.unpack(\">II\", f.read(8))\n",
    "    label = np.fromfile(f, dtype=np.dtype(np.uint8).newbyteorder('>'))\n",
    "\n",
    "\n",
    "data = np.array(data, dtype=float) / 255\n",
    "label = np.array(label, dtype=int)\n",
    "\n",
    "print('array shapes:', data.shape, label.shape)\n",
    "\n",
    "# number of samples\n",
    "n_smpl = data.shape[0]\n",
    "n_px = data.shape[1]\n",
    "n_feat = n_px**2\n",
    "\n",
    "# format in scikit learn format\n",
    "X = data.reshape([n_smpl,-1])\n",
    "y = label\n",
    "\n",
    "# number of classes\n",
    "n_cl = np.unique(y).size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c89d8cd9-0c27-4050-b77d-2fc18b0a21ea",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# plot example\n",
    "plt.figure()\n",
    "plt.imshow(data[0], cmap='binary')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45051d6f",
   "metadata": {},
   "source": [
    "We want to compare the accuracy and reverse image of different classifiers:\n",
    "- MLR is the multinomial logistic regression\n",
    "- 1NN is the 1-nearest-neighbor\n",
    "- RF is a random forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ffeb39f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# classifiers\n",
    "\n",
    "clf = LogisticRegression(penalty='l2', C=0.1, max_iter=1000)\n",
    "\n",
    "# clf = KNeighborsClassifier(n_neighbors=3, algorithm='brute', metric='minkowski')\n",
    "\n",
    "# clf = RandomForestClassifier(n_estimators=10, criterion='gini')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ed9e3c",
   "metadata": {},
   "source": [
    "Usual cross-validation scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08581b70",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 80% for training and 20% for testing\n",
    "cvs = StratifiedShuffleSplit(n_splits=20, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc99137f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# dataframe to store result accuracies\n",
    "df_acc = pd.DataFrame()\n",
    "\n",
    "# repeat classification\n",
    "for train_ind, test_ind in cvs.split(X, y):\n",
    "    # fit to train data\n",
    "    clf.fit(X[train_ind,:], y[train_ind])\n",
    "    # test accuracy\n",
    "    d = {'type': ['test'],\n",
    "         'score': [clf.score(X[test_ind,:], y[test_ind])]}\n",
    "    df_acc = pd.concat((df_acc, pd.DataFrame(data=d)), ignore_index=True)\n",
    "\n",
    "    # shuffling labels and fit again for evaluation of chance level\n",
    "    train_ind_rand = np.random.permutation(train_ind)\n",
    "    clf.fit(X[train_ind,:], y[train_ind_rand])\n",
    "    d = {'type': ['shuf'],\n",
    "         'score': [clf.score(X[test_ind,:], y[test_ind])]}\n",
    "    df_acc = pd.concat((df_acc, pd.DataFrame(data=d)), ignore_index=True)\n",
    "        \n",
    "# print table of results\n",
    "print(df_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77102885",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# test versus shuffle accuracy\n",
    "sb.violinplot(data=df_acc, x='type', y='score', order=['shuf','test'],\n",
    "              density_norm='width', palette=['gray','orange'])\n",
    "plt.plot([-1,2], [1.0/n_cl]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2368ba7-584b-47c4-8b9c-ea6836daf72b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# input PCA space\n",
    "n_comp = 10\n",
    "pca = PCA(n_components=n_comp)\n",
    "pca.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c21965d9-ded9-452e-a5c4-c48dc71ea48a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# train classifier\n",
    "clf.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8e740a4-fb72-4e0f-9a9c-9c55ea916727",
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_bub = 2.0\n",
    "prec_bub = np.eye(2) / sig_bub**2\n",
    "det_prec_bub_2pi = 1.0 / 2 / np.pi / sig_bub**2\n",
    "\n",
    "def gaussian_2d(x, mean, prec, det_prec_2pi):\n",
    "    # assert means.shape==(2,) and covs.shape==(2,2)\n",
    "    return np.exp(-np.einsum('kli,ij,klj->kl', x-mean, prec, x-mean) / 2.0) * det_prec_2pi\n",
    "\n",
    "def gen_bubbles(n_bub=10):\n",
    "    x, y = np.meshgrid(np.arange(n_px), np.arange(n_px))\n",
    "    pos = np.dstack((x, y))\n",
    "    means = np.random.randint(n_px, size=[n_bub,2])\n",
    "    xy = np.zeros([n_px,n_px])\n",
    "    for i_bub in range(n_bub):\n",
    "        xy += gaussian_2d(pos, means[i_bub,:], prec_bub, det_prec_bub_2pi)\n",
    "    return xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c679b749-4e86-4979-b19f-8a51ea27b521",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(gen_bubbles().reshape([n_px,n_px]))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "183e7d57-9a6b-46a4-93fa-f16b83635b09",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# number of noisy samples\n",
    "n_ns = 10000\n",
    "\n",
    "# noisy inputs\n",
    "method = 'raw' # raw noise, pca noise, bubble, alone or combined with stimulus\n",
    "if method=='raw':\n",
    "    # noise in original image space\n",
    "    X_ns = np.random.normal(loc=0.0, scale=1.0, size=[n_ns,n_feat])\n",
    "elif method=='pca':\n",
    "    # noise in PCA space\n",
    "    X_ns = pca.inverse_transform(np.random.normal(loc=0.0, scale=2.0, size=[n_ns,n_comp]))\n",
    "elif method=='pca+stim':\n",
    "    # additive noise in PCA space on top of stimulus\n",
    "    X_ns = X[np.random.randint(n_smpl, size=[n_ns]),:]\n",
    "    X_ns += pca.inverse_transform(np.random.normal(loc=0.0, scale=0.5, size=[n_ns,n_comp]))\n",
    "elif method=='bub':\n",
    "    # bubbles are Gaussian randomly distributed in original space (2D image)\n",
    "    bub_masks = np.zeros([n_ns,n_px,n_px])\n",
    "    for i_ns in range(n_ns):\n",
    "        bub_masks[i_ns,:,:] = gen_bubbles()\n",
    "    bub_masks = bub_masks.reshape([n_ns,n_px**2])\n",
    "    X_ns = bub_masks * 50.0\n",
    "elif method=='bub*stim':\n",
    "    # mulitplicative noise bubble\n",
    "    X_ns = X[np.random.randint(n_smpl, size=[n_ns]),:]\n",
    "    bub_masks = np.zeros([n_ns,n_px,n_px])\n",
    "    for i_ns in range(n_ns):\n",
    "        bub_masks[i_ns,:,:] = gen_bubbles()\n",
    "    bub_masks = bub_masks.reshape([n_ns,n_px**2])\n",
    "    X_ns *= bub_masks * 100.0\n",
    "else:\n",
    "    raise ValueError('unknown method')\n",
    "\n",
    "# prediction from trained classifier\n",
    "y_pred = clf.predict(X_ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3fa6c2a-536e-4d69-9fb3-9b582c63fc65",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# check correct prediction for digit 0\n",
    "i_cl = 3\n",
    "y_correct = np.array(y_pred==i_cl, dtype=int)\n",
    "print('ratio correct predicitons:', y_correct.sum()/n_ns)\n",
    "\n",
    "# correlation between prediction for digit 0 and inputs\n",
    "X_y_1 = np.einsum('si,s->i', X_ns, y_correct) # 1st moment\n",
    "X_y_2 = np.einsum('si,s->i', X_ns**2, y_correct) # 2nd moment\n",
    "X_corr = ( X_y_2 - X_y_1**2 ) / n_ns # covariance\n",
    "X_corr /= np.std(X_ns, axis=0) * np.std(y_correct)\n",
    "\n",
    "#  plot \"reverse image\"\n",
    "data_inv = -X_corr.reshape([n_px,n_px])\n",
    "plt.figure()\n",
    "plt.imshow(data_inv, cmap='binary')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}