hue_map.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "animal = 'AE2'\n",
    "unit = '001'\n",
    "trial = '021'\n",
    "\n",
    "file_path = 'D:\\\\2P_data\\\\F\\\\%s\\\\u%s_%s\\\\1\\\\F_%s_u%s_%s_plane1_proc.mat' % (animal, unit, trial, animal, unit, trial)\n",
    "analyzer_path = 'F:\\\\NHP\\\\AE2\\\\AnalyzerFiles\\\\%s\\\\%s_u%s_%s.analyzer' % (animal, animal, unit, trial)\n",
    "plot = 'y'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.io as sio\n",
    "import functions\n",
    "\n",
    "F = sio.loadmat(file_path, squeeze_me=True, struct_as_record=False)['dat']\n",
    "params = functions.analyzer_params(analyzer_path)\n",
    "trial_num, stim_time = functions.analyzer_pg_conds(analyzer_path)\n",
    "trial_num['direction'] = trial_num.ori\n",
    "trial_num.ori[(trial_num.ori >= 180) & (trial_num.ori != 256)] = trial_num.ori[(trial_num.ori >= 180) &\n",
    "                                                                               (trial_num.ori != 256)] - 180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "mean_image = F.mimg[:, :, 1]\n",
    "color_mask = np.zeros(mean_image.shape)\n",
    "pref_cells = pd.DataFrame(columns=['expt', 'x', 'y', 'color_pref'])\n",
    "osi_df = pd.DataFrame(columns=['expt', 'x', 'y', 'r_pref', 'r_orth', 'osi', 'color_pref'])\n",
    "color_hex = ['#af1600', '#8a4600', '#5a5d01', '#2a6600', '#006a00', '#006931', '#006464', '#0058b6', '#002DFF', '#6a2ade', '#97209b', '#aa1c50']\n",
    "\n",
    "for i, cell in enumerate(F.stat):\n",
    "    if cell.iscell == 1:\n",
    "        scans = int(F.Fcell.shape[1] / trial_num.shape[0])\n",
    "        t = np.linspace(stim_time[0] * -1, stim_time.sum() - stim_time[1], scans)\n",
    "        cell_F = pd.DataFrame(np.reshape(F.Fcell[i, :], (int(F.Fcell.shape[1] / scans), scans)))\n",
    "        cell_F.insert(0, 'Orientation', trial_num.ori)\n",
    "        cell_F.insert(1, 'Direction', trial_num.direction)\n",
    "        cell_F.insert(2, 'Spatial Frequency', trial_num.s_freq)\n",
    "        cell_F.insert(3, 'Color', trial_num.colormod)\n",
    "        cell_F.insert(4, 'Baseline',\n",
    "                      np.reshape(F.Fcell[i, :], (int(F.Fcell.shape[1] / scans), scans))\n",
    "                      [:, 0:int(np.floor(scans * stim_time[0] / stim_time.sum())) - 1].mean(axis=1))\n",
    "        \n",
    "        start_sample = int(np.ceil(scans * (stim_time[0] + 0.5) / stim_time.sum()))\n",
    "        end_sample = int(scans - np.ceil(scans * (stim_time[1]) / stim_time.sum()))\n",
    "        \n",
    "        cell_F.iloc[:, -scans:] = cell_F.iloc[:, -scans:].subtract(cell_F['Baseline'], axis=0).divide(cell_F['Baseline'],\n",
    "                                                                                                  axis=0)\n",
    "        \n",
    "        cell_F_stim = pd.DataFrame(cell_F.iloc[:, -(scans-start_sample):-(scans-end_sample)].mean(axis=1), columns=['Mean'])\n",
    "        cell_F_stim.insert(0, 'Orientation', trial_num.ori)\n",
    "        cell_F_stim.insert(1, 'Direction', trial_num.direction)\n",
    "        cell_F_stim.insert(2, 'Spatial Frequency', trial_num.s_freq)\n",
    "        cell_F_stim.insert(3, 'Color', trial_num.colormod)\n",
    "        \n",
    "        max_response = cell_F_stim.groupby('Color').mean()['Mean'].max()\n",
    "        \n",
    "        # Calculate OSI and save into DataFrame\n",
    "        x_med, y_med = cell.med\n",
    "        counts_ori = cell_F_stim.groupby(['Orientation', 'Spatial Frequency']).mean()\n",
    "        counts_sf = cell_F_stim.groupby(['Spatial Frequency', 'Orientation']).mean()\n",
    "        sem_ori = cell_F_stim.groupby(['Orientation', 'Spatial Frequency']).sem()\n",
    "        sem_sf = cell_F_stim.groupby(['Spatial Frequency', 'Orientation']).sem()\n",
    "        \n",
    "        s_freq_max = counts_ori.unstack(level=1).mean(axis=0).idxmax()\n",
    "        orientation_max = counts_ori.unstack(level=1).mean(axis=1).idxmax()\n",
    "        ori_pref = counts_sf['Mean'][s_freq_max[1]].idxmax()\n",
    "        r_pref = counts_sf['Mean'][s_freq_max[1]].loc[ori_pref]\n",
    "        if ori_pref < 90:\n",
    "            ori_orth = ori_pref + 90\n",
    "        else:\n",
    "            ori_orth = ori_pref - 90\n",
    "        r_orth = counts_sf['Mean'][s_freq_max[1]].loc[ori_orth]\n",
    "        osi = (r_pref - r_orth) / r_pref\n",
    "        \n",
    "        osi_df.loc[i] = [analyzer_path[-17:-9], x_med, y_med, r_pref, r_orth, osi, 0]\n",
    "        \n",
    "        colors = np.unique(trial_num.colormod.values)\n",
    "        list_F_stim = [cell_F_stim[cell_F_stim['Color'] == color]['Mean'] for color in colors]\n",
    "        F_value, p = sp.stats.f_oneway(*list_F_stim)\n",
    "        \n",
    "        if p < 0.001:\n",
    "            max_response = cell_F_stim.groupby('Color').mean()['Mean'].max()\n",
    "            if max_response > 0.04:\n",
    "                max_color = cell_F_stim.groupby('Color').mean()['Mean'].idxmax()\n",
    "                if max_color < 16:\n",
    "                    ortho_color = max_color + 6\n",
    "                else:\n",
    "                    ortho_color = max_color - 6\n",
    "                ortho_response = cell_F_stim.groupby('Color').mean()['Mean'][ortho_color]\n",
    "                color_si = (max_response - ortho_response) / max_response\n",
    "                if color_si > 0.7:\n",
    "                    color_mask[(cell.ypix - 1), (cell.xpix - 1)] = max_color\n",
    "                    if plot == 'y':\n",
    "                        plt.ioff()\n",
    "                        fig, ax = plt.subplots()\n",
    "                        ax.scatter(cell_F_stim.groupby('Color').mean()['Mean'].index.values, cell_F_stim.groupby('Color').mean()['Mean'].values, color=color_hex, s=100, zorder=2)\n",
    "                        ax.fill_between(cell_F_stim.groupby('Color').mean()['Mean'].index.values, \n",
    "                                        cell_F_stim.groupby('Color').mean()['Mean'] - cell_F_stim.groupby('Color').sem()['Mean'],\n",
    "                                        cell_F_stim.groupby('Color').mean()['Mean'] + cell_F_stim.groupby('Color').sem()['Mean'],\n",
    "                                        alpha=0.5, color='#808080', zorder=1)\n",
    "                        ax.xaxis.set_ticks([])\n",
    "                        ax.set_ylabel(r'Response ($\\Delta F/F$)')\n",
    "                        plt.savefig('D:/nhp_cells/cell%d.pdf' % i, format='pdf')\n",
    "                        plt.close()\n",
    "                    pref_cells.loc[i] = [analyzer_path[-17:-9], x_med, y_med, max_color]                    \n",
    "                    osi_df.loc[i, 'color_pref'] = 1\n",
    "\n",
    "pref_cells.to_csv('D:/nhp_cells/%s_color_pref_df.csv' % analyzer_path[-17:-9])\n",
    "osi_df.to_csv('D:/nhp_cells/color_osi/%s_color_osi_df.csv' % analyzer_path[-17:-9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "from matplotlib.colors import ColorConverter\n",
    "# colors = ['#BA0000', '#C03100', '#D3ED00', '#3AF000', '#09F300', '#0AF125', '#0FEBAD', '#032EB3', '#000091', '#23009B', '#1200A2', '#B50020']\n",
    "colors = ['#af1600', '#8a4600', '#5a5d01', '#2a6600', '#006a00', '#006931', '#006464', '#0058b6', '#002DFF', '#6a2ade', '#97209b', '#aa1c50']\n",
    "cmap = ListedColormap(colors, 'indexed')\n",
    "cmap.set_under(ColorConverter.to_rgba('white', alpha=0))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 7))\n",
    "ax.imshow(mean_image, cmap='gray', alpha=0.8)\n",
    "im = ax.imshow(color_mask, cmap=cmap, vmin=10)\n",
    "cb = fig.colorbar(im)\n",
    "cb.set_ticks([])\n",
    "cb.set_label('Preferred Color')\n",
    "plt.axis('off')\n",
    "plt.savefig('D:/nhp_cells/%s_color_map.pdf' % analyzer_path[-17:-9], format='pdf', bbox_inches='tight')\n",
    "np.save('D:/nhp_cells/%s_color_map.npy' % analyzer_path[-17:-9], color_mask)\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ4AAAD7CAYAAACbmXq7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF9NJREFUeJzt3X2wHWV9B/DvN5AQIAFKEuNoIwmQ0iYiaU0rgiN0OogB\npDeOsYVgalucxAl9w0KtooBjHesbFCcQICjlpUHTkusLEdpaUEFbe4NMMlABhYlBLCZUEMJLEvLr\nH7t7s2fv7p7n2d1zdvec72fmDPecu7tnc7nne3/Ps88+D80MIiI+JtV9AiLSPgoOEfGm4BARbwoO\nEfGm4BARbwoOEfGm4CiB5Kkkn6j7PLoheTTJ5+s+DxkcCg4AJM8lOUbyeZI/I/kNkm/p8zk8GL7/\n8yRfIflS7PmHyhzbzB4zs2lVnWscyVtIGskzE69/Pnz9vF68r9Rr6IOD5IUArgTwCQCzAbwOwBoA\nZ/f4fQ+IPzezhWY2LfyAfwfABdFzM/tEL8+lAo8AWBE9ITkZwLsAPFbbGUlPDXVwkDwcwMcArDaz\n281sl5ntMbOvm9nF4TYHkbyS5JPh40qSB2Uc7zdI3kPymbCCODv2vRtJXkNyE8ldAH7X81w/TvLG\n2PNjSVrs+b0kLyf5XZLPkbyT5JG+24bf/2OSPyG5k+SHSD5B8tSc0xsFcGr48wSAMwGMAdgRO+Z8\nkneT/L/wuDfHtkf4Hn9D8n9I/oLkDVk/Z6nfUAcHgDcDmApgY842HwZwIoBFAE4A8DsALkluFP6V\n/RqAfwXwKgB/BuBWksfFNjsXwN8BmA7g3grOP+lcAH+EoHI6FMCFvtuSPB7AVQD+EMBrAcwC8Oou\n7/sigDsAvDt8vgLATYltCODj4bEWADgawEcS2ywHcBqA+QAWAvjbLu8rNRn24JgBYKeZ7c3ZZjmA\nj5nZz81sB4DLAbwnZbsTAUwD8Ekz221m/wHg6wDOiW3zFTO7z8z2mdlLFf0b4m4ws0fN7AUAGxCE\nne+2ywCMmtl3zexlpIRkhpsArAgrl5MAfDX+TTN7xMy+Gf5sfg7gCgCnJI5xlZk9YWY7ETQdz4E0\n0oF1n0DNngYwk+SBOeHxGgDbYs+3ha+lbbfdzPYltn1t7Pn2Mifr4H9jX7+AIMh8t30NYudpZrtI\n/sLhvb8F4FcRVAlfMbOXSY5/k+SrEVQyJyOouCYh1pQJxX8+WT9naYBhrzi+B+BlACM52zwJ4KjY\n89eFr6VtN4fkpMS2P409L3Mr8i4Ah8Sed2s+FPUzBAEAACB5KIBf6baTBbdZ3wrgA5jYTAGAv0fw\nsz7ezA4D8F4EzZe4ObGvs37O0gBDHRxm9iyAjwJYQ3KE5CEkJ5NcQvJT4WbrAVxCchbJmeH2t6Qc\n7r8Q/OW+ODzGqQDeAeC2ik73AQCnkJxD8ggAH6zouEkbAIyQPJHkFASdx66uAHCamd2X8r3pCMLv\nWZJzAPx1yjYXkHwtyRkIKpcveZ679MlQBwcAmNlnEXQMXoKgdN4O4AIEVwqAoENvDMAWAFsB3B++\nljzObgRBsQTATgBXA1hhZj+s6FTvRNCJuxXA95HoQ6iKmW0B8FcIAuRJBM25pxFUC932fdrMvpnx\n7UsRdCw/i+Dc/yVlm/UA/h3AjwE8jKCfQxqImshH8pA8DMAzAI4ys5710YQjcM8zs3t69R5SnaGv\nOGQikmeHzbZpAD4L4P5ehoa0j4JD0ixF0Ex5AsBc6LKoJKipIiLeVHGIiDcFh4h48xo5OnPmTJs7\nd27Hazu27szdZ9bxM71PSmSQbd68eaeZzSpzjEVTZ9tz+3Y7bfvYnmfuMrO3l3m/JK/gmDt3LsbG\nxsafr523ruug4FVj5xc6MZFBRXJb963yPbdvNz45+1Snbd/9xGjlf71LNVVWPZ4dCqsePz/3+yLS\nXqX7ONICQoEhMtgq7xxVaIgMvspuq1dgiAyPnlyO5RuCh4gMpkqDIxkYCg+RwVRpcNiWia+p+hAZ\nPIX6OIoEwdp56zK/p/4RkXbpy5yj1zyn0JDy4n+w0qpb6Z++3Kvy/ukKBykurbmrJnC9CgVHkbRP\nqyxUbUg3CodmqrypYluyS0oFRXVGZ++fIHzkKc2pIv1VODjyAkLtz96Kh4ZIHbyaKpsfmhgW0UOk\nF/J+t/R7V59Sl2P1P64ew9Y00e9Z82gGMBHxNuxrx4q00l4ejqcPPNNx69Hum3gqXHGofBQZXl7B\n8cYF6gwVkYrvVVGgiAyHSvs4dLVFeiHtBsnkbQz6neuvyq+q6H+gVM1lxLGGpvdXoYojOWo0eq0M\nNX/ElW6arF+pIedVyftrwTcoPGR/1bF23rrxaRqSAaLflf5pxTgOzcMgQP5kUNJftY8c9W2bah6G\n4RVVHe+ffn5qc0V/VPqn9oojrb9EqsdzvtHx3NYvqelMyln1+PlYBVWhdau94ihKYeMvHhbJIGkb\n3Zldr9orDmB/1RH9EriEgn5h3MQDIv51WysOaYbGVByaCKi/2l5xSL0aExxJCo9q2Pol448khYcU\n1YimiiuFSTlqnkhVGh0cCgqRZmp0cMjg2TBnaerry7Zv7POZSBmN7eOQ4ZIVKNJMqjikMTbMWarK\nw9GLk6Zj6yG/V9v7q+IQEW+qOKSRNFit2VRxSF8t277RuznCc76hMScNo4pDarFs+8bxDtG0IImq\njLQh86pA6qfgkNq4VB5pASL1U3BIK8QDRP0f9VMfh7RKMijU/1EPBYeIeFNTRVpHzZP6qeIQEW+q\nOCRTVt+B/uKLKg5JldfhqM5IUcUhjcX3vS/1dbv++j6fiSSp4pBGygoNaQYFh6Rqcj+GQqV+aqpI\nJlu/JLU/IwqVtMl32jKfRrQEh6anLEbBIbmyKo+sGbvaNhmPFqouRsEhlet1eJTpHE1b7KuN4bGD\nB+HqA+bV9v7q45BG6ueVk7aFRhMoOKSQvIqiqmrDrr9+QoBkBQpX0+2YiRUDFRrFqKkihUUBsXAl\ncdmmEVx2xigevNYqf59u1UcUGvHwsDXZ56GwKE8Vh1TisjNGAQQhUre80JBqqOIoIHlFoU1XEaqW\nFhQLV7InlUcWBUX/qeLwlHYZcsOcpVpQSIaKgkNKefBam1Bd9LPakHqoqeIpPjt38vVhprAYLgqO\nAqKQyJveX7pL9o8ofNpDTZUSiiwuJIGsTtV+v7/r+A/ppOCQRqnjci5XUwHiScEhfVf3WI/o/ZPn\nofBwpz4OcTIsHcJReGhsSD5VHNJV1WNU8jpB+9FBGl1Cfijnz6aqj3wKDulqECuLhSuJBXsnvm5r\nbPwh2dRUkVp1qzDS/vInP9S+zYt438aCvW5VTnIej2G/UU7BIU6yBr4V5fRhzWguqBlRDskDANxk\nZsuLHkPBIc7i4TEozZci1UYjTH0WWOi4vs3Wzqdm9grJo0hOMbPdRd5ewSFemhwYXE2n5opGqAIA\nHgNwH8mvAtgVvWhmn3PZWcEhjWVrLLOPo1szpsrOzar7M445EG+s9oiF/Dh8TAIw3XdnBYc0WlYA\n5IVHlaLlIapYZ2Z0dnP6ZszscgAgOS18/rzP/gqOARPvwOzWrPDZtkmyqop4kLg2W3LfJ7amTHJ9\nGd8gaVJoAADJ1wO4GcCR4fOdAFaY2YMu+2scx4BIm0wo7ypI2rZtm4yoLfeYNC00QtcBuNDMjjKz\nowB8AIDz1PKqOFou+rBP3j3a8fqeKSN1nI60x6Fmdnf0xMzuIXmo686qOFosKzSi19Jedz1mm8Sb\nLvHmSVSR1F2VjDzVyKs4j5H8CMm54eMSBFdanCg4QmvnrcPaeevqPg0vvZoPpOnhkdZ30YtwsPVL\nKlt8u4Hh8ScAZgG4PXzMCl9zQjP3f9DixYttbGzM9wRbIxkcqx4/v6YzcZdXdQDZv7A+4RAPp+Ur\ngr81t960L7Xt3sAPSOOQ3Gxmi0sd48j5htOvctv4tjNKv1+S+jhi2hAUSdFozj1TRryaJmWHkC9f\nMQnLUl6PwqRogLjcmyLFkfwagMwfqJmd7XIcBUfNqlyjJdkh2u1YyblTs2QtIr3hdGLZXcU+1EXW\nXomGfg/7DWYlfaaKgyg4atD0PoSkKDSiZkoZ8TtTk+GR108Rv1+kjavLN4WZfSv6muQUAL8WPn3Y\nzPa4Hkedo+Jc5dx6077O/QpWG3ELVxILV+Zf+eBqTlgsWsoheSqARwGsAXA1gEdIvtV1fwXHgPJt\n8mRtn/l6BaFRhEKjMp8F8DYzO8XM3grgdABXuO6spsoAyAsJn2HlLmFz6037wmP6XcWJPHitpU5W\nnBx/kaYfoZF3C/2AhdZkM3s4emJmj5Cc7LqzgqPlkpdEow9u1hq3kXhIFOmgjXfERivVA8CHK1hw\nuq6rKI2cd6N3xkiuA3BL+Pw8AM5jLRQcNXC9mtFN2jiK0dl0vhSa9f7xqyijs9l1+Pplm0Y6wiOq\nKHwCJNnx+sDBE9ek7SWX0BiwTtn3A1gN4M/D598GcI3rzgqOGhW99NqPm6aqeI+sS67J16q4WlOW\nbekeHo0KjQOeAabd7r0byVkAZpnZQwA+Fz5AciGAwwDscDlO/f/HxIvLB9qlknGtdsrcLOdSMeSF\nRnTFpe4FnAbM5wHMTHn9SAD/4HoQBUfLZDVDRp6y8YdLJdNtG9/AuGxT5/auzYzkJV4gaKYk9SM8\nbEt2VdGoaqOcY83s28kXzew7AJx7edRUkUx5Ha1J8T6OMn0bWfrZ3zFAIZEmb5pAXVUZZC6dn2Wq\njnhQ5IVGfP+0+1ZcJCuO5SsmYdGLzN2mSln9GgMcHj8ieYaZbYq/SHIJPG6rV3AMMJ8rLBGf5Q96\nMfVgL0MiKa8zdMCuoMT9JYA7SL4bwObwtcUA3gzgLNeDKDgGRPxDXGQCn+QxokuyTx94Zuq2M/be\nUeg90iT7L7xvfssaMNaD8SBtv9HOzB4leTyAcwG8Pnz5WwBWmtlLrsdRcAyAvOZE1gAx1+POAFLD\n4+kDzxwPjzonOu56j0uF4TEoN9qZ2csAvljmGLqq0nJpoZF3RSTvcm5WAMzYewdm7L2jY76SKkOj\nDQskpTVr+IahG206ThXHAKqiqRLxvfkt6dOn7w+qizJujEs2VYrM1ZElr+pwGfTl9B4trj6KUnAM\niSJXYtKCxHWWtHhgJF9LBkj8xreqq49uTZXoAz+MH/4yFBwtl7YQtO89LMmxFNGVjV70XWRVHXU1\nV5L9FsBgBwjJrcifOtCpBlMfxwBIznaeDIhoRGmatAFYVdw7khUQaZVIGd0qCtdZxbq9nhcmLQua\nswC8A8Cd4WN5+NgUPpwoOAbUyFOGPVNGsGfKSO5Vl6xxE7248eyiuywzUJomKzzSHm1iZtvMbBuA\n08zsYjPbGj4+COBtrsdRcAwJ3yUel91lpe+QjQdF0wJjWK+GxJDkybEnJ8EjD9THMaBc5troF9fQ\n8BmJGg+1jQCWvsvvnKq6otJifwrgCyQPD58/A48FmVRxDCjXcBidPXGJg/jzuhZMzquQGrqIc6uY\n2WYzOwHACQBOMLNFZna/6/6qOFoq/uGJd3zmLYiUFiYjT1lqeADAo4s6q4Wsjs2yzRDfmdCic47b\n+M+d3493imaN5cirOtrWd+GL5GwAnwDwGjNbQnIBgDeb2Q0u+ys4WiotLLrJuoEt+UF8dNH+73W7\nCuITGvHlI11lNa2if39WgLoONW/tOA7uAiZv7r5dthsRDDv/cPj8EQBfAuAUHGqqtNTobI4/fKX9\nhc+7ZJvH5fLq8hWTOq7SxL/uVm30qz+mVaFRjZlm9mUA+wDAzPYCeMV1Z1UcA6roRMjxaqOXlq+Y\nNHGBpwKTOGuR68J2kZyBcDAYyRMBPOu6s4JjQJVdVLpOva4y4v0fQ7yg9YUAvgrgGJL3AZgFwPna\nlIJjgPU6PPL6N+IfznORf0dudI4+kwgVkTeKdJiQnARgKoBTABwHgPBcO1bB0VJpnYNp0sr/Kj6Y\nrqGRZuTu3+84p7Ln4zLniEJjPzPbR3KNmf0mgAeLHINm7qXa4sWLbWzMebEnaamsDs+5D3XO8xH/\nwKftc/GxwX8/9aPs94of0zdA8kIzHh55odHPpkp0B/BjtwMv7rBSScZXHWJ416+7bXzNDzab2eKO\n/cnPAPgegNvNJwRCqjiktKygyQuMpEFvnjRwbZiVCPo59pJ8CUFzxczsMJedFRwyQVozJG992Yvu\nskJ3vUbVRtWh0ZQrLQ0MCwDBTSoAFprZT4oeQ+M4xFsVH/R4E6VIB25WM6UpoQE0d0rEsGlSarZp\nVRzSIesDuazHH0jfjtIiI2eTet2/kVVxvLQTpYZ8VuR+kr9tZv9dZGcFhzjJuwcG6D70PN6USXay\nlpVXZUThUMfYjXjFEQ+RqTPxxr6cQL43AVhOchuAXdjfx+F0z7CCQ3oiWQXM7/xuqcWsi6h7oFc8\nRHgdm1BxnF5mZwWH1GLy7s7wqHNtlmFkZttIvgXAfDP7IslZAKa57q/gGHBZHY9ZH9S0W9Z9+Ow7\nefcoNoRNmNGUe1ekd0heimDpx+MQ3CU7GcAtAE7O2y+iqypSGZ/QeHRR8FgU65/wmee0zN3BAgBY\nCuBsBP0bMLMnkb+SfQdVHAMu634VnykEu13iHO/4DO+snf+A1yli0VOGByqocpIT+AD192002G4z\nM5LR3bGH+uys4JAJyo6FSN6anxYk8x+YOMOYq6wqI22UaNXrxw6QL5O8FsARJN+HYL7R6113VnAM\nqao6I8usk1J0ysFkP8x40GUMLx/M8HgR4A+89yJ5kJm9bGafIXkagF8i6Of4qJn9m+txFBxDoExI\nuIzf6BYeLhVM8hjdQiUeHtEqdckxG4MXFpX4HoDfInmzmb0HgHNYxCk4JFP8r3peJ+RFFYwqLVJ9\nJAMpOseNKd+TcVNIngvgJJLvTH7TzG53OYiCQybwvVIRX5c2bd86PsR5a+UOuVUIlnw8AsFSkHEG\nQMEh/dPEy6IKj4nM7F4A95Icc10KIY3GcchAa2KgNcRtJC8heR0AkJxP8izXnRUc0qEXH7R+fXjT\nqgtVHJm+AGA3gJPC5z8F8HHXndVUkQ7dhpyndUhmdVJm7dNLCgpnx5jZH5A8BwDM7IVwgh8nCg5x\n5voXXR/eVthN8mDsX1flGAAvu+6s4JAJ4ldIFAID61IAdwKYQ/JWBDe3vdd1Z/VxSCaX0Fg7b13H\n1/Hn0kxhk+SHAN6JICzWA1hsZve4HkMVhxQWhUQyLNbOW4dVj5/f9/PxHX06rMKb2zaZ2fEoOPeo\ngkMKi8Khrioj3gmbtubtp0+nwiNbqTlH1VSRUrJCo5dhkjYPh++t/II3AfhPkj8muYXkVpJbXHdW\nxSGlRFVHUztSVXVk0pyjIsD+OT4kG8mpCO5XORbAVgA3mNle3+MoOKSwZHOh2y34/aDw6OofAewB\n8B0ASwAsAPAXvgdRcAyZKj7cTb7/IxkaaqZMsCC8mgKSNwD4fpGDKDhkYA1yaEw14GjHBsZDnU/3\nRF+Y2V6PUeYdFBziLet+ln41UZrYCdsiJ5D8Zfg1ARwcPtdq9ZKuyiaGPrztZGYHVHEcjeMYIvqw\nS1VUcQwZhYdUQRWHiHhTcIiINwWHiHhTcIiIN3WOSlfxRaurWjpS2k0Vh4h4U8UhmeKVhkicKg7J\npGaJZFHFIbkUHpJGFYeIeFNwiIg3BYeIeFNwiIg3dY6KtNAh+4BFL7rNr/IQqr8jWhWHiHhTxTGk\nuLrzr5Wt0Twd4k4Vh4h4U8UxpIa1wkgOo9cAt2JUcchQ2zBnqe7JKUDBIUNDAVEdNVVaos51TESS\nVHG0QNZ6KE1eirGJ1J9RHQWHDBWFRzUUHA2nqqJ68fBYtn2jwqQA9XE0XNSPkddcUV+HP4VFOao4\nWkLhIE2i4Gg5BYrUQU2VFlFISFOo4hARbwoOEfGm4BARbwoOEfGmzlGRFjp4L7DIsbP8n3rw/qo4\nRMSbgkNEvCk4RMSbgkNEvCk4RMSbgkNEvCk4RMSbgkNEvCk4RMSbgkNEvCk4RMSbgkNEvCk4RMSb\ngkNEvCk4RMSbgkNEvCk4RMSbZgATaaGpLwDzH6jv/VVxiIg3BYeIeFNwiIg3BYeIeFNwiIg3BYeI\neFNwiIg3BYeIeFNwiIg3BYeIeFNwiIg3BYeIeFNwiIg3BYeIeFNwiIg3BYeIeKOZuW9M7gCwrXen\nIzIUjjKzWWUOQPJOADMdN99pZm8v834T3t8nOEREADVVRKQABYeIeFNwiIg3BYeIeFNwiIg3BYeI\neFNwiIg3BYeIeFNwiIi3/weLqOR/yqookQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x20f0878a4a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import glob\n",
    "color_map_files = glob.glob('D:\\\\nhp_cells\\\\*color_map.npy')\n",
    "fig, ax = plt.subplots()\n",
    "for file in color_map_files:\n",
    "    color_map = np.load(file)\n",
    "    ax.imshow(color_map, cmap=cmap, vmin=10, vmax=21)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "cb = fig.colorbar(im)\n",
    "cb.set_ticks([])\n",
    "cb.set_label('Preferred Color')\n",
    "plt.title('Color Tuning Map')\n",
    "plt.savefig('D:/nhp_cells/color_map_all.pdf', format='pdf', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_pref_files = glob.glob('H:\\\\nhp_cells\\\\*color_pref_df.csv')\n",
    "list_ = []\n",
    "for file in color_pref_files:\n",
    "    df = pd.read_csv(file, index_col=0, header=0)\n",
    "    list_.append(df)\n",
    "color_pref_all = pd.concat(list_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_pref_rb = color_pref_all[(color_pref_all.color_pref == 17) | (color_pref_all.color_pref == 10)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import Normalize\n",
    "from sklearn import svm\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import StratifiedShuffleSplit, GridSearchCV\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(-3, 3, 500),\n",
    "                     np.linspace(-3, 3, 500))\n",
    "X = color_pref_rb[['x', 'y']].as_matrix()\n",
    "scaler = StandardScaler()\n",
    "X = scaler.fit_transform(X)\n",
    "Y = color_pref_rb[['color_pref']].as_matrix().ravel()\n",
    "\n",
    "clf = svm.NuSVC()\n",
    "clf.fit(X, Y)\n",
    "\n",
    "Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
    "Z = Z.reshape(xx.shape)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(Z, interpolation='nearest', extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
    "          origin='lower', cmap='jet_r')\n",
    "contours = ax.contour(xx, yy, Z, levels=[-0.5, .5], linewidths=1, colors='k')\n",
    "scat = ax.scatter(X[:, 0], X[:, 1], s=30, c=Y, cmap=cmap, edgecolors='k', vmin=10, vmax=21)\n",
    "fig.colorbar(scat)\n",
    "plt.savefig('H:/nhp_cells/color_map_all_fit.pdf', format='pdf', bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_osi_files = glob.glob('H:\\\\nhp_cells\\\\color_osi\\\\*color_osi_df.csv')\n",
    "list_ = []\n",
    "for file in color_osi_files:\n",
    "    df = pd.read_csv(file, index_col=0, header=0)\n",
    "    list_.append(df)\n",
    "color_osi = pd.concat(list_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFO9JREFUeJzt3X+MXeV95/H3d8egEifL7dbdCtmwtlakHXdSVukU+ses\n6lmrxZRNSZSuljFplDCR5TSMslqtZUtTBbFhEjsoUrcqieVlLJPdMPxT4pLYMVmJuaVWwgazAYI9\nAlmQGEMkL0lFMg5aGPPdP+ZiXU/mx73jc+feO+f9kkaa+5xn7vnq6DEfzjnPeU5kJpKk8vpn7S5A\nktReBoEklZxBIEklZxBIUskZBJJUcgaBJJWcQSBJJWcQdJGI+BcR8Y2IOB8RP46I7e2uSWqViLgr\nIk5ExP+LiEPtrmc1W9PuAtSU+4G3gN8C/g1wJCKezcyT7S1LaonXgHuBm4Gr2lzLqhY+WdwdImIt\n8E9AX2a+WGv7GvBaZu5pa3FSC0XEvcCGzPxEu2tZrbw01D3eD8y8GwI1zwK/26Z6JK0SBkH3eC/w\n8zltPwfe14ZaJK0iBkH3mAb++Zy2q4FftKEWSauIQdA9XgTWRMT1dW03AN4olnRZDIIukZnngUeA\n/xoRayNiAPgz4H+0tzKpNSJiTUT8GtAD9ETEr0WEMx1bwCDoLn/J7DS6c8BDwKedOqpV7K+AN4E9\nwMdqv/9VWytapZw+Kkkl5xmBJJWcQSBJJWcQSFLJGQSSVHIGgSSVXNvm5K5bty43btzYrt13tfPn\nz7N27dp2l9F1nn766dcz8zfbtX/H/PI43pev0THftiDYuHEjJ06caNfuu1q1WmXLli3tLqPrRMSP\n27l/x/zyON6Xr9Ex76UhSSo5g0CSSs4gkKSSMwgkqeQMAkkqOYOgi0xMTNDX18fWrVvp6+tjYmKi\n3SVJWgVc27tLTExMMDo6yvj4OBcuXKCnp4fh4WEAhoaG2lydpG7mGUGXGBsbY3x8nMHBQdasWcPg\n4CDj4+OMjY21uzRJXW7JIIiIgxFxLiKeX6TPloh4JiJORsQ/FFuiAKamphgYGLikbWBggKmpqTZV\nJGm1aOSM4BCwbaGNEVEBvgL8WWb+LvAfiilN9Xp7ezl+/PglbcePH6e3t7dNFUlaLZYMgsx8AvjZ\nIl22A49k5pla/3MF1aY6o6OjDA8PMzk5yczMDJOTkwwPDzM6Otru0iR1uSJuFr8fuCIiqsD7gP+W\nmV+br2NE7AB2AZVKpUK1Wi1g9+VwzTXXcMcdd3DnnXdy5swZrrvuOj72sY9xzTXXeBw7mGP+8k1P\nT3vcWqyhdxZHxEbgW5nZN8+2vwX6ga3Mvlj9e8CtmfniYt/Z39+fLsC1PC7CtTwR8XRm9rdr/475\n5XG8L1+jY76IM4KzwE8z8zxwPiKeAG4AFg0CSVJnKGL66N8DAxGxJiLeA9wEOJVFkrrEkmcEETEB\nbAHWRcRZ4G7gCoDM3J+ZUxFxDHgOeAd4IDMXnGoqSeosSwZBZi752Gpm3gfcV0hFkqQV5ZPFXcS1\nhiS1gmsNdQnXGpLUKp4RdAnXGpLUKgZBl3CtIUmtYhB0CdcaktQqBkGXcK0hSa3izeIu8e4N4ZGR\nEaampujt7WVsbMwbxZIum0HQRYaGhhgaGnLtFUmF8tKQJJWcQSBJJWcQSFLJGQSSVHIGgSSVnEEg\nSSVnEEhSyS0ZBBFxMCLORcSiL5uJiD+IiJmI+PPiypMktVojZwSHgG2LdYiIHmAf8J0CapIkraAl\ngyAznwB+tkS3EeDvgHNFFCVJWjmXfY8gItYDHwG+evnlSJJWWhFrDf01sDsz34mIRTtGxA5gF1Cp\nVCpUq9UCdl8+09PTHrsu4Zi/fI731ovMXLpTxEbgW5nZN8+2l4F3E2Ad8EtgR2YeXuw7+/v788SJ\nE83WK3DRuWWKiKczs79d+3fML4/jffkaHfOXfUaQmZvqdnqI2cBYNAQkSZ1jySCIiAlgC7AuIs4C\ndwNXAGTm/pZWJ0lquSWDIDMbfvNJZn7isqqRJK04nyyWpJIzCCSp5AwCSSo5g0CSSs4g6GARMe/P\n4ODggtuk1WJiYoK+vj62bt1KX18fExMT7S5p1SriyWK1yEIP+23cc4Qf7b11hauRVs7ExASjo6OM\nj49z4cIFenp6GB4eBmBoqOGJjGqQZwSSOs7Y2Bjj4+MMDg6yZs0aBgcHGR8fZ2xsrN2lrUoGgaSO\nMzU1xcDAwCVtAwMDTE1Ntami1c0gkNRxent7OX78+CVtx48fp7e3t00VrW4GgaSOMzo6yvDwMJOT\nk8zMzDA5Ocnw8DCjo6PtLm1V8maxpI7z7g3hkZERpqam6O3tZWxszBvFLWIQSOpIQ0NDDA0NuQz1\nCvDSkCSVnEEgSSVnEEjqSD5ZvHIaeTHNQeDfA+cWeFXlHcBuZl9X+Qvg05n5bNGFSioPnyxeWY2c\nERwCti2y/WXgjzLzA8DngQMF1CWpxMbGxti+fTsjIyPcfPPNjIyMsH37dp8sbpFG3lD2RO3l9Qtt\n/27dxyeBDZdflqQyO3XqFOfOnWPt2rVkJufPn+fAgQO8/vrr7S5tVSp6+ugw8O2Cv1NSyfT09DAz\nM8PBgwcvXhr66Ec/Sk9PT7tLW5UKC4KIGGQ2CAYW6bMD2AVUKpUK1Wq1qN2XjseuOzjml2dmZgaA\nZ555hk2bNvHyyy9fbPcYFi8WWur4kk6zl4a+Nd/N4tr23wO+AdySmS82suP+/v48ceJE45XqIpeh\nXp6IeDoz+9u1f8d84yKCPXv28M1vfvPik8Uf+tCH2Lt374LLs+tXNTrmL/uMICKuAx4B/qLREJCk\nxWzYsIEHH3yQr3/96xcvDd1xxx1s2OAtyFZoZProBLAFWBcRZ4G7gSsAMnM/8DngN4Cv1N6QNdPO\n/+uS1P2+9KUv8dnPfpY777yTM2fOcN111zEzM8OXv/zldpe2KjUya2jRSbuZ+SngU4VVJKn03n1W\n4N3pomvXruULX/iCzxC0iIvOSepILjq3clxiQpJKziCQpJIzCCSp5AwCSSo5g0CSSs4gkKSSMwgk\nqeQMAkkqOYNAkkrOIJCkkjMIJKnkDAJJKjmDQJJKziCQpJJbMggi4mBEnIuI5xfYHhHxNxFxOiKe\ni4gPFl+mJKlVGjkjOARsW2T7LcD1tZ8dwFcvvyxJ0kpZMggy8wngZ4t0uQ34Ws56EqhExDVFFShJ\naq0i7hGsB16p+3y21iZJ6gIr+qrKiNgB7AIqlUqFarW6krtfVTx23cExf/mmp6c9bi1WRBC8Clxb\n93lDre1XZOYB4ABAf39/+h7SZTp2xHe4dgnH/OXzncWtV8SloUeBj9dmD/0h8EZm/qSA75UkrYBG\npo9OAN8DfjsizkbEcETsjIidtS5HgZeA08B/B/6yZdVKKo2JiQn6+vrYunUrfX19TExMtLukVWvJ\nS0OZObTE9gQ+U1hFkkpvYmKC0dFRxsfHuXDhAj09PQwPDwMwNLTof5K0DD5ZLKnjjI2NMT4+zuDg\nIGvWrGFwcJDx8XHGxsbaXdqqZBBI6jhTU1MMDAxc0jYwMMDU1FSbKlrdVnT6qCQ1ore3l3vuuYfD\nhw8zNTVFb28vH/7wh+nt7W13aauSQSCp4wwODrJv3z727dvH5s2bOXXqFLt372bnzp1L/7GaZhBI\n6jiTk5Ps3r2bgwcPXjwj2L17N4cPH253aauSQSCp40xNTfGDH/yAe++99+IDZW+//TZf/OIX213a\nqmQQSOoIEXHJ5yuvvHLJfrOz13W5nDUkqSNk5sWfhx56iE2bNvH4449z3X85zOOPP86mTZt46KGH\nLumnYnhGIKnjvPvQ2MjICGdOTTHy7V7GxsZ8mKxFDAJJHWloaIihoSE27jnC83tvbXc5q5qXhiSp\n5AwCSSo5g0CSSs4gkKSSMwgkqeQaCoKI2BYRL0TE6YjYM8/2qyPimxHxbEScjIhPFl+qJKkVGnlD\nWQ9wP3ALsBkYiojNc7p9BjiVmTcAW4AvR8T8jwVKkjpKI2cENwKnM/OlzHwLeBi4bU6fBN4Xs89+\nvxf4GTBTaKWSpJZoJAjWA6/UfT5ba6v3t0Av8BrwQ+CzmflOIRVKklqqqCeLbwaeAf4d8K+B/xUR\n/5iZP6/vFBE7gF1ApVKpUK1WC9p9+XjsuoNjvhget9ZqJAheBa6t+7yh1lbvk8De2ovsT0fEy8Dv\nAN+v75SZB4ADAP39/blly5Zlll1yx47gsesOjvkCON5brpFLQ08B10fEptoN4NuBR+f0OQNsBYiI\n3wJ+G3ipyEIlSa2x5BlBZs5ExF3AY0APcDAzT0bEztr2/cDngUMR8UMggN2Z+XoL65YkFaShewSZ\neRQ4Oqdtf93vrwF/UmxpkqSV4JPFklRyBoEklZxBIEklZxBIUskZBJJUcgaBJJWcQSBJJWcQSFLJ\nGQSSVHIGgSSVnEEgSSVnEEhSyRX1Yhpdhhvu+Q5vvPl2U3+zcc+RhvtefdUVPHu3awJKmp9B0AHe\nePNtfrT31ob7V6vVpl7U0UxoSCofLw1JUsk1FAQRsS0iXoiI0xGxZ4E+WyLimYg4GRH/UGyZkqRW\nWfLSUET0APcDfwycBZ6KiEcz81RdnwrwFWBbZp6JiH/ZqoIlScVq5IzgRuB0Zr6UmW8BDwO3zemz\nHXgkM88AZOa5YsuUJLVKI0GwHnil7vPZWlu99wO/HhHViHg6Ij5eVIGSpNYqatbQGuD3ga3AVcD3\nIuLJzHyxvlNE7AB2AZVKpUK1Wi1o992vmWMxPT3d9LHzWLeHY74YHrfWaiQIXgWurfu8odZW7yzw\n08w8D5yPiCeAG4BLgiAzDwAHAPr7+7OZKZCr2rEjTU0HbXb6aLPfr+I45gvg+G25Ri4NPQVcHxGb\nIuJK4Hbg0Tl9/h4YiIg1EfEe4CZgqthSJUmtsOQZQWbORMRdwGNAD3AwM09GxM7a9v2ZORURx4Dn\ngHeABzLz+VYWLkkqRkP3CDLzKHB0Ttv+OZ/vA+4rrjRJ0krwyWJJKjmDQJJKziCQpJIzCCSp5AwC\nSSo5g0CSSs4gkKSSMwgkqeQMAkkqOYNAkkrOIJCkkjMIJKnkinoxjSQ15IZ7vsMbb77d1N9s3HOk\n4b5XX3UFz979J82WVWoGgaQV9cabb/Ojvbc23L/ZFzE1Exqa5aUhSSq5hoIgIrZFxAsRcToi9izS\n7w8iYiYi/ry4EiVJrbRkEERED3A/cAuwGRiKiM0L9NsHfKfoIiVJrdPIGcGNwOnMfCkz3wIeBm6b\np98I8HfAuQLrkyS1WCNBsB54pe7z2VrbRRGxHvgI8NXiSpMkrYSiZg39NbA7M9+JiAU7RcQOYBdQ\nqVQqVKvVgnbf/Zo5FtPT000fO491ezjm5+d47yyNBMGrwLV1nzfU2ur1Aw/XQmAd8KcRMZOZh+s7\nZeYB4ABAf39/NjMlbFU7dqSp6XHNTqdr9vtVHMf8PBzvHaeRIHgKuD4iNjEbALcD2+s7ZOamd3+P\niEPAt+aGgCSpMy0ZBJk5ExF3AY8BPcDBzDwZETtr2/e3uEZJUgs1dI8gM48CR+e0zRsAmfmJyy9L\nkrRSfLJYkkrOIJCkkjMIJKnkXH20A7yvdw8feHDBJZzm92Az3w/Q+GqPksrFIOgAv5ja67K8ktrG\nS0OSVHIGgSSVnEEgSSVnEEhSyXmzWNKKcpZc5zEIJK0oZ8l1Hi8NSVLJGQSSVHIGgSSVnEEgSSXX\nUBBExLaIeCEiTkfEr9zuj4g7IuK5iPhhRHw3Im4ovlRJUissGQQR0QPcD9wCbAaGImLznG4vA3+U\nmR8APk/tHa2SpM7XyBnBjcDpzHwpM98CHgZuq++Qmd/NzH+qfXyS2RfcS5K6QCPPEawHXqn7fBa4\naZH+w8C3L6eoMmp67vOxxvtffdUVTVYjqUwKfaAsIgaZDYKBBbbvAHYBlUqlQrVaLXL3XevQtrVN\n9f/EsfNN/43Huj0c8/Nr5jhMT083fdw8zs1pJAheBa6t+7yh1naJiPg94AHglsz86XxflJkHqN0/\n6O/vz2aeFlSdY0eaetJS7eOYn0eT47fZJ4v999G8Ru4RPAVcHxGbIuJK4Hbg0foOEXEd8AjwF5n5\nYvFlSpJaZckzgsyciYi7gMeAHuBgZp6MiJ217fuBzwG/AXwlIgBmMrO/dWVLkorS0D2CzDwKHJ3T\ntr/u908Bnyq2NEnSSnD1UUkrzllyncUgkLSimlmCGmZDo9m/UXNca0iSSs4gkKSSMwgkqeQMAkkq\nOYNAkkrOIJCkkjMIJKnkDAJJKjmDQJJKziCQpJIzCCSp5AwCSSo5g0CSSs7VRztY7SU/82/bN397\nZraoGkmrVUNnBBGxLSJeiIjTEbFnnu0REX9T2/5cRHyw+FLLJzPn/ZmcnFxwmyQ1a8kgiIge4H7g\nFmAzMBQRm+d0uwW4vvazA/hqwXVKklqkkTOCG4HTmflSZr4FPAzcNqfPbcDXctaTQCUirim4VklS\nCzRyj2A98Erd57PATQ30WQ/8pL5TROwAdgGVSqVCtVpttl4B09PTHrsu4Zhv3ODg4ILbFronNjk5\n2aJqymVFbxZn5gHgAEB/f39u2bJlJXe/alSrVTx23cEx37iF7nE53luvkUtDrwLX1n3eUGtrto8k\nqQM1EgRPAddHxKaIuBK4HXh0Tp9HgY/XZg/9IfBGZv5k7hdJkjrPkpeGMnMmIu4CHgN6gIOZeTIi\ndta27weOAn8KnAZ+CXyydSVLkorU0D2CzDzK7H/s69v21/2ewGeKLU2StBJcYkKSSs4gkKSSMwgk\nqeQMAkkquWjXQmUR8X+BH7dl591vHfB6u4voQv8qM3+zXTt3zC+b4335GhrzbQsCLV9EnMjM/nbX\nIa0Ex3vreWlIkkrOIJCkkjMIutOBdhcgrSDHe4t5j0CSSs4zAkkqOYNgFYiIB+Z5fai0ajnmi+Wl\nIUkqOc8IOlhE/OeIeL72858iYm1EHImIZ2tt/7HWrxoRzrNW13PMt8eKvqpSjYuI32f2vQ43AQH8\nb2bfB/FaZt5a63N1+yqUiuWYbx/PCDrXAPCNzDyfmdPAI8DbwB9HxL6I+LeZ+UZ7S5QK5ZhvE4Og\n+3wQ+CFwb0R8rt3FSCvAMd9iBkHn+kfgwxHxnohYC3wEeBr4ZWb+T+A+Zv+BSKuFY75NvEfQoTLz\n/0TEIeD7taYHgPcC34+Id5g9Zf50m8qTCueYbx+nj0pSyXlpSJJKziCQpJIzCCSp5AwCSSo5g0CS\nSs4gkKSSMwgkqeQMAkkquf8Pd5DPpTquIZgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x19582031390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grouped_osi = color_osi[color_osi.r_pref > 0.02][['osi', 'color_pref']].groupby(['color_pref'])\n",
    "grouped_osi.boxplot(column=['osi'])\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.20052198772076169"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_stat, p = stats.ttest_ind(color_osi[(color_osi.color_pref == 1) & (color_osi.r_pref > 0.03)]['osi'],\n",
    "                            color_osi[(color_osi.color_pref == 0) & (color_osi.r_pref > 0.03)]['osi'])\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'color_osi' is not defined",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-81-ee96b498940f>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mcolor_osi\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'color_osi' is not defined"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    ""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}