notebook_local_not_precomputed.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "\n",
    "import ipywidgets as ipy\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from PIL import Image\n",
    "from scipy.stats import zscore\n",
    "from skimage import io\n",
    "from sklearn.cluster import DBSCAN\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.manifold import TSNE\n",
    "from sklearn.metrics import accuracy_score, make_scorer\n",
    "from sklearn.model_selection import cross_val_predict, KFold, LeaveOneOut, GridSearchCV\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.preprocessing import minmax_scale\n",
    "from sklearn.svm import SVC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. read single image, and display it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### read image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_img = io.imread('/path/to/timepoint_12.tiff')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### check data type and shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(single_img))\n",
    "print(single_img.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### plot single image interactively (it is an image stack that has 40 layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@ipy.interact(layer=ipy.widgets.IntSlider(value=26, min=0, max=39))\n",
    "def f(layer):\n",
    "    fig, ax = plt.subplots(figsize=(4, 8)) # here you can change figure size\n",
    "    ax.matshow(single_img[layer,:,:], cmap='hot')\n",
    "\n",
    "# in case you want to play with colormaps:\n",
    "# https://matplotlib.org/examples/color/colormaps_reference.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### in case you had troubles with the previous widget, you can simply use following cell (not interactive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(4, 8)) # here you can change figure size\n",
    "ax.matshow(single_img[28,:,:], cmap='hot') # integer number is layer index (here 28)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need data in two different forms\n",
    "1. `images_dict`. It it dictionary (Python data structure)\n",
    "2. `images_matrix`. 2-dimensional numpy array\n",
    "\n",
    "Simply run the following cell to prepare data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_template = '/path/to/data/'\n",
    "images_dict = dict()\n",
    "images_matrix = np.zeros([100, 40*256*180])\n",
    "i = 0\n",
    "for img_name in sorted(os.listdir(path_template)):\n",
    "    img = io.imread(os.path.join(path_template, img_name))\n",
    "    images_dict[img_name[:-5]] = img\n",
    "    images_matrix[i] = np.reshape(img, [-1], order='C')\n",
    "    i = i+1\n",
    "\n",
    "print(images_matrix.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Work towards PCA analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to perform PCA analysis, images need to be preprocessed in the certain way. This part of the notebook will guide you through this process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. subtract background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### interactive widget\n",
    "Here you can interactively explore which threshold value is most appropriate.\n",
    "1. threshold slider -> to pick different values of the threshold\n",
    "2. img_layer slider -> to browse different layer in the images stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@ipy.interact(threshold=ipy.widgets.IntSlider(value=90, min=90, max=170))\n",
    "def f(threshold):\n",
    "    timepoint = np.copy(images_dict['timepoint_12'])\n",
    "    timepoint[timepoint <= threshold] = 0\n",
    "\n",
    "    @ipy.interact(img_layer=ipy.widgets.IntSlider(value=26, min=0, max=39))\n",
    "    def f(img_layer):\n",
    "        fig, ax = plt.subplots(figsize=(5,10))\n",
    "        ax.matshow(timepoint[img_layer,:,:], cmap='hot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### non-interactive version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 120 # manipulate this\n",
    "timepoint = np.copy(images_dict['timepoint_12'])\n",
    "timepoint[timepoint <= threshold] = 0\n",
    "fig, ax = plt.subplots(figsize=(5,10))\n",
    "ax.matshow(timepoint[26,:,:], cmap='hot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 compare before-after"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@ipy.interact(img_layer=ipy.widgets.IntSlider(value=26, min=0, max=39))\n",
    "def f(img_layer):\n",
    "    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(10,10))\n",
    "    ax1.matshow(images_dict['timepoint_12'][img_layer,:,:], cmap='hot')\n",
    "    ax1.set_title('before', size='26')\n",
    "    ax2.matshow(timepoint[img_layer,:,:], cmap='hot')\n",
    "    ax2.set_title('after', size='26')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(10,10))\n",
    "ax1.matshow(images_dict['timepoint_12'][26,:,:], cmap='hot')\n",
    "ax1.set_title('before', size='26')\n",
    "\n",
    "ax2.matshow(timepoint[26,:,:], cmap='hot')\n",
    "ax2.set_title('after', size='26')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. remove pixels with zeros"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 get indices with nonzero pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Procedure is as follows:\n",
    "1. vectorize image stack `timepoint`.\n",
    "2. get only nonzero positions in the vector\n",
    "\n",
    "As a result the `inx` stores indices of nonzero values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timepoint_vector = np.reshape(timepoint, [-1], order='C')\n",
    "timepoint_vector.shape\n",
    "\n",
    "inx = np.nonzero(timepoint_vector)\n",
    "print(inx[0])\n",
    "print(len(inx[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 get only foreground pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images_matrix_fg = np.zeros([100, len(inx[0])]) # fg = foreground\n",
    "for i in range(images_matrix.shape[0]):\n",
    "    images_matrix_fg[i,:] = images_matrix[i,inx[0]]\n",
    "\n",
    "print('images matrix has shape:            {}'.format(images_matrix.shape))\n",
    "print('images matrix foreground has shape: {}'.format(images_matrix_fg.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. normalize input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images_matrix_fg_norm = zscore(images_matrix_fg, axis=1)\n",
    "images_matrix_fg_norm = np.nan_to_num(images_matrix_fg_norm)\n",
    "\n",
    "#im2d_n = zscore(im2d_fg, axis=0) # im2d_n = image 2d normalized (double check the axis (0,1,None))\n",
    "#im2d_n = np.nan_to_num(im2d_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 check normalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('max value:   {}'.format(images_matrix_fg_norm.max()))\n",
    "print('min value:   {}'.format(images_matrix_fg_norm.min()))\n",
    "print('mean values: {}'.format(images_matrix_fg_norm.mean(axis=1)))\n",
    "print('std values:  {}'.format(images_matrix_fg_norm.std(axis=1)))\n",
    "#im2d_n.std(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. clip outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images_matrix_fg_norm_clipped = np.clip(images_matrix_fg_norm, a_min=-3.0, a_max=3.0)\n",
    "print('min: {}'.format(images_matrix_fg_norm_clipped.min()))\n",
    "print('max: {}'.format(images_matrix_fg_norm_clipped.max()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. calculate PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 50\n",
    "pca = PCA(n_components=n_components)\n",
    "pca_space = pca.fit_transform(images_matrix_fg_norm_clipped)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1 check output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('output data type: ' + str(type(pca_space)))\n",
    "print('output shape:     ' + str(pca_space.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.1 prepare principal components for visualization\n",
    "Reshape back to the original size (256, 180, 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_space_bg = np.zeros([40*256*180, n_components]) # pca_space_bg = PCA space with background (zeros)\n",
    "for i in range(pca_space_bg.shape[1]):\n",
    "    pca_space_bg[inx, i]= pca.components_[i,:] # Note that here we make use of the previously computed inx\n",
    "pca_space_bg = pca_space_bg.reshape([40, 256, 180, -1])\n",
    "print(pca_space_bg.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. visualize principal components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### interactive widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@ipy.interact(principal_component=ipy.widgets.IntSlider(value=0, min=0, max=59))\n",
    "def f(principal_component):\n",
    "    @ipy.interact(img_layer=ipy.widgets.IntSlider(value=26, min=0, max=39))\n",
    "    def f(img_layer):\n",
    "        fig, ax = plt.subplots(figsize=(7,7))\n",
    "        cax = ax.matshow(pca_space_bg[img_layer,:,:, principal_component], cmap='PiYG')\n",
    "        fig.colorbar(cax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interactive wigdet: display only max projections (easier to visually inspect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_abs = np.abs(pca_space_bg)\n",
    "@ipy.interact(principal_component=ipy.widgets.IntSlider(value=0, min=0, max=59))\n",
    "def f(principal_component):\n",
    "    fig, ax = plt.subplots(figsize=(7,7))\n",
    "    ax.matshow(pca_abs.max(axis=0)[:,:,principal_component], cmap='hot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### non-interactive: display only max projections (easier to visually inspect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_abs = np.abs(pca_space_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "ax.matshow(pca_abs.max(axis=0)[:,:,0], cmap='hot') # change int here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# three plot on one figure\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(12,12))\n",
    "ax1.matshow(pca_abs.max(axis=0)[:,:,0], cmap='hot') # change int here (range = 0, n_components)\n",
    "ax2.matshow(pca_abs.max(axis=0)[:,:,1], cmap='hot') # change int here (range = 0, n_components)\n",
    "ax3.matshow(pca_abs.max(axis=0)[:,:,2], cmap='hot') # change int here (range = 0, n_components)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. visualize variance explained ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(10,5), facecolor='w')\n",
    "ax1.plot(pca.explained_variance_ratio_, color='black', lw=2)\n",
    "ax1.grid(axis='both')\n",
    "ax1.set_ylabel('proportion of variance explained', fontsize=10)\n",
    "ax1.set_xlabel('principal component', fontsize=10)\n",
    "\n",
    "ax2.plot(np.cumsum(pca.explained_variance_ratio_), color='black', lw=2)\n",
    "ax2.grid(axis='both')\n",
    "ax2.set_ylabel('cumulative proportion of variance explained', fontsize=10)\n",
    "ax2.set_xlabel('number of components', fontsize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### end of the PCA analysis\n",
    "\n",
    "If this task was easy for you, you may want to try different dimensionality reduction methods. Sklearn offers coherent API for many methods.\n",
    "Please refer to sklearn.decomposition module: http://scikit-learn.org/stable/modules/classes.html#module-sklearn.decomposition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# calculate t-SNE embeddings on the PCA space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.manifold import TSNE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### prepare PCA space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 99\n",
    "pca = PCA(n_components=n_components)\n",
    "pca_space = pca.fit_transform(images_matrix_fg_norm_clipped)\n",
    "print('pca space dimensions: {}'.format(pca_space.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### explore t-SNE embeddings as a function of the perplexity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@ipy.interact(perplexity=ipy.widgets.IntSlider(value=29, min=0, max=50))\n",
    "def f(perplexity):\n",
    "    @ipy.interact(n_iter=ipy.widgets.IntSlider(value=1000, min=500, max=5000))\n",
    "    def g(n_iter):\n",
    "        tsne = TSNE(n_components=2, n_iter=n_iter, perplexity=perplexity)\n",
    "        tsne_output = tsne.fit_transform(pca_space)\n",
    "        print('t-SNE output has shape: {}'.format(tsne_output.shape))\n",
    "        fig, ax = plt.subplots(figsize=(7,7))\n",
    "        ax.scatter(tsne_output[:,0], tsne_output[:,1])\n",
    "        ax.grid(axis='both')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### non-interactive version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perplexity=29 # manipulate this\n",
    "n_iter=1000\n",
    "tsne = TSNE(n_components=2, n_iter=n_iter,  perplexity=perplexity)\n",
    "tsne_output = tsne.fit_transform(pca_space)\n",
    "print('t-SNE output has shape: {}'.format(tsne_output.shape))\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "ax.scatter(tsne_output[:,0], tsne_output[:,1])\n",
    "ax.grid(axis='both')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# cluster t-SNE embeddings using DBSCAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import DBSCAN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### explore clustering results as a function of two parameters: eps and nim_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsne = TSNE(n_components=2, n_iter=1000,  perplexity=29)\n",
    "tsne_output = tsne.fit_transform(pca_space)\n",
    "@ipy.interact(eps=ipy.widgets.FloatSlider(value=1.2, min=0, max=10))\n",
    "def f(eps):\n",
    "    @ipy.interact(min_samples=ipy.widgets.IntSlider(value=5, min=3, max=50))\n",
    "    def g(min_samples):\n",
    "        dbscan = DBSCAN(eps=eps, min_samples=min_samples)\n",
    "        dbscan_output = dbscan.fit_predict(tsne_output)\n",
    "        plt.scatter(x=tsne_output[:,0], y=tsne_output[:,1], c=dbscan_output)\n",
    "        plt.grid()\n",
    "        print(dbscan_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### non-interactive version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps=1.9 # manipulate this\n",
    "min_samples=5 # manipulate this\n",
    "\n",
    "tsne = TSNE(n_components=2, n_iter=1000,  perplexity=29)\n",
    "tsne_output = tsne.fit_transform(pca_space)\n",
    "\n",
    "dbscan = DBSCAN(eps=eps, min_samples=min_samples)\n",
    "dbscan_output = dbscan.fit_predict(tsne_output)\n",
    "plt.scatter(x=tsne_output[:,0], y=tsne_output[:,1], c=dbscan_output)\n",
    "plt.grid()\n",
    "print(dbscan_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Developmental state classification\n",
    "Early vs. late development of the zebrafish\n",
    "labeling scheme:\n",
    "* early = 0\n",
    "* late  = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 15\n",
    "pca = PCA(n_components=n_components)\n",
    "pca_space_late = pca.fit_transform(images_matrix_fg_norm_clipped)\n",
    "print('pca space dimensions: {}'.format(pca_space.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_space_early = np.load('/path/to/pca_space_early.npy', allow_pickle=False)\n",
    "print(type(pca_space_early))\n",
    "print(pca_space_early.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### prepare input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.concatenate((np.zeros([100]).astype(np.int), np.ones([100]).astype(np.int)), axis=0)\n",
    "X = np.concatenate((pca_space_early[:,0:5], pca_space_late[:,0:5]), axis=0)\n",
    "\n",
    "print(X.shape)\n",
    "print(y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### prepare experiment with SVC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### linear kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 3 # manipulate this\n",
    "cv_scheme = KFold(n_splits=10, shuffle=True, random_state=120715)\n",
    "model = SVC(C=C, kernel='linear')\n",
    "y_pred = cross_val_predict(model, X, y, cv=cv_scheme, n_jobs=1, verbose=0, method='predict')\n",
    "print(accuracy_score(y, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `rbf` kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 3 # manipulate this\n",
    "cv_scheme = KFold(n_splits=10, shuffle=True, random_state=120715)\n",
    "model = SVC(C=3, kernel='rbf')\n",
    "y_pred = cross_val_predict(model, X, y, cv=cv_scheme, n_jobs=1, verbose=0, method='predict')\n",
    "print(accuracy_score(y, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualize timepoints with highest probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timepoint_early = np.load('/path/to/timepoint_68_early.npy', allow_pickle=False)\n",
    "timepoint_late = np.load('/path/to/timepoint_63_late.npy', allow_pickle=False)\n",
    "print(type(timepoint_early))\n",
    "print(type(timepoint_late))\n",
    "print(timepoint_early.shape)\n",
    "print(timepoint_late.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(10,10))\n",
    "ax1.matshow(timepoint_early, cmap='hot')\n",
    "ax1.set_title('early', size='26')\n",
    "ax2.matshow(timepoint_late, cmap='hot')\n",
    "ax2.set_title('late', size='26')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### run grid search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SVC()\n",
    "parameters = {'kernel':('linear', 'rbf'), 'C': [0.01,1,3,7]}\n",
    "scoring=make_scorer(accuracy_score)\n",
    "cv_scheme = KFold(n_splits=10, shuffle=True, random_state=120715)\n",
    "gs = GridSearchCV(model, parameters, scoring=scoring, cv=cv_scheme)\n",
    "gs.fit(X,y)\n",
    "cv_results_df = pd.DataFrame(gs.cv_results_)\n",
    "cv_results_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### prepare experiment with Gaussian Naive Bayes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_scheme = KFold(n_splits=10, shuffle=True, random_state=120715)\n",
    "model = GaussianNB()\n",
    "y_pred = cross_val_predict(model, X, y, cv=cv_scheme, n_jobs=1, verbose=0, method='predict')\n",
    "print(accuracy_score(y, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final remarks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Final chart with some insights on one visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_chart = pca_abs.max(axis=0)[:,:,[2,9,12]]\n",
    "final_chart = final_chart / final_chart.max()\n",
    "final_chart = final_chart*255\n",
    "img = Image.fromarray(final_chart.astype(np.uint8))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "ax.imshow(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# end"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}