{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Make sure retinatools is in the python path\n",
    "import sys\n",
    "sys.path.append('../') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab as pl\n",
    "from retinatools.library import filter_video_PC, get_PC_response, filter_video_MC, get_MC_response, filter_video_BO, get_BO_response, compare_to_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: all reliabilities rounded to 2 significant digits, and rounded up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.rcParams['figure.figsize'] = [15, 4]\n",
    "def make_plots(cell_name, RpcF, RpfS, filename, dt, shift, column, retinal_reliability): \n",
    "    \n",
    "    # Calculate the R^2 values / variaces explained\n",
    "    rF = compare_to_file(filename, RpcF, dt, column, shift)    \n",
    "    rS = compare_to_file(filename, RpcS, dt, column, shift)\n",
    "    spatial_factors = rS**2\n",
    "    temporal_factors = rF**2 - rS**2\n",
    "    missing_in_model = retinal_reliability - rF**2\n",
    "    \n",
    "    #Get the data for the plot\n",
    "    cell_activity = np.loadtxt(filename)\n",
    "    time_cell = np.arange(0, len(cell_activity))*dt + 5000 - shift*6.66 \n",
    "    time_model = np.arange(0, len(RpcS)*dt, dt)\n",
    "\n",
    "    #Make the plot\n",
    "    pl.figure()\n",
    "    pl.subplot(1,2,1)\n",
    "    pl.plot(time_model, RpcS/np.std(RpcS), label=\"Model without surround\")\n",
    "    pl.plot(time_model, RpcF/np.std(RpcF), label = \"Full model\")\n",
    "    pl.plot(time_cell, cell_activity[:,column]/np.std(cell_activity[:,column]), label=\"Cell data\")\n",
    "    pl.xlim([10000,11500])\n",
    "    pl.ylim([0,8])\n",
    "    pl.xlabel(\"Time [ms]\")\n",
    "    pl.ylabel(\"Activity [Z-scored]\")\n",
    "    pl.legend();\n",
    "    pl.title(\"r_model=\" + str(rF))\n",
    "    pl.subplot(1,2,2)\n",
    "    labels = 'Center-contribution', '\"Surround\"-Contribution', 'Missing in Model', 'Individual variations'\n",
    "    sizes = [spatial_factors, temporal_factors, missing_in_model, 1-retinal_reliability]\n",
    "    explode = (0.1, 0.1, 0.1, 0.1)\n",
    "    pl.pie(sizes, explode=explode, labels=labels, autopct='%1.1f%%', shadow=True, startangle=-60)\n",
    "    pl.axis('equal')\n",
    "    pl.title(\"Cell: \" +  cell_name);\n",
    "\n",
    "    #Save figure + data\n",
    "    pl.savefig(\"./\"+ cell_name + \".pdf\")\n",
    "    \n",
    "    # Select data for cell and model, and save it\n",
    "    t0 = 5000\n",
    "    time_cell = np.arange(0, len(cell_activity))*dt + 5000 - shift*6.66\n",
    "    cell = cell_activity[(time_cell>t0) * (time_cell<64980), column]\n",
    "    time = np.arange(0, len(RpcF)*dt, dt) \n",
    "    model = RpcF[(time>=t0) * (time<64980)]\n",
    "    time_file = time[(time>=t0) * (time<64980)]\n",
    "    model_data_save = np.transpose([time_file, cell/np.std(cell), model/np.std(model)])\n",
    "    np.savetxt(\"./\"+ cell_name + '_time_cell_fit.csv', model_data_save, fmt='%1.5f', delimiter=', ', newline='\\n', header='time [ms], cell, model', footer='', comments = '', encoding=None)\n",
    "\n",
    "    pl.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PC cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cell in [\"110Ron\", \"130Gon\"]:\n",
    "    \n",
    "    if cell == \"110Ron\":\n",
    "        fitparameters = {'xoff': 2.9, 'yoff': -1.8, 'rcenterM': 20, 'rcenterL': 5.0, \\\n",
    "                         'lum_factor': 1200, 'tau1': 6.9, 'tau2': 60, 'k1': 0.017, 'tau3': 400,\\\n",
    "                         'alpha': 0.68, 'k2': 4, 'o1': 0.07, 'dt': 1000/150.}\n",
    "        filename = './data/File110Ron.txt'\n",
    "        celltype = \"Ron\"\n",
    "        fitparameters_s = fitparameters.copy()\n",
    "        fitparameters_s['rcenterM'] = fitparameters_s['rcenterL']\n",
    "        retinal_reliability = 0.91**2 # From Barry for 110Ron\n",
    "\n",
    "    elif cell == \"130Gon\":\n",
    "        fitparameters = {'xoff': 2.9, 'yoff': -2.2, 'rcenterM': 4.9, 'rcenterL': 32, \\\n",
    "                         'lum_factor': 290, 'tau1': 7.8, 'tau2': 100, 'k1': 0.018, 'tau3': 2000,\\\n",
    "                         'alpha': 0.71, 'k2': 3.0, 'o1': 0.06, 'dt': 1000/150.}\n",
    "        filename = './data/File130Gn.txt'\n",
    "        celltype = \"Gon\"\n",
    "        fitparameters_s = fitparameters.copy()\n",
    "        fitparameters_s['rcenterL'] = fitparameters_s['rcenterM']\n",
    "        retinal_reliability = 0.80**2 # From Barry for 130Gon\n",
    "\n",
    "\n",
    "    # full model\n",
    "    lum_signalM, lum_signalL = filter_video_PC('../stimuli/1x10_256.mpg', 9750, fitparameters)\n",
    "    RpcF = get_PC_response(lum_signalM, lum_signalL, fitparameters, celltype)\n",
    "\n",
    "    # no-surround model\n",
    "    lum_signalM, lum_signalL = filter_video_PC('../stimuli/1x10_256.mpg', 9750, fitparameters_s)\n",
    "    RpcS = get_PC_response(lum_signalM, lum_signalL, fitparameters_s, celltype)\n",
    "    \n",
    "    make_plots(cell, RpcF, RpcS, filename, dt = fitparameters['dt'], shift = 1, column = 0, retinal_reliability = retinal_reliability)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MC cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cell in [\"299Mon\", \"078Moff\"]:\n",
    "    \n",
    "    if cell == \"299Mon\":\n",
    "        fitparameters = {'xoff': -2.5, 'yoff': -5.5, 'rcenter': 4.0, 'rsurr': 13.0, 'ksurr': 0.45, 'surround_delay': 3, \\\n",
    "                         'lum_factor': 60, 'tau1': 7.3, 'tau2': 38, 'k1': 0.019, 'c1': 0.016, 'tauplus': 110, \\\n",
    "                         'tauA': 20, 'k2': 50, 'dt': 1000/150.}\n",
    "        filename = './data/File299Mon.txt'\n",
    "        celltype = 'on-cell'\n",
    "        fitparameters_s = fitparameters.copy()\n",
    "        fitparameters_s['rsurr'] = fitparameters_s['rcenter']\n",
    "        retinal_reliability = 0.92**2 # From 299 MC on cell\n",
    "\n",
    "    elif cell == \"078Moff\":\n",
    "        fitparameters = {'xoff': -8.4, 'yoff': 2.5, 'rcenter': 4.5, 'rsurr': 20, 'ksurr': 0.34, 'surround_delay': 8, \\\n",
    "                         'lum_factor': 1000, 'tau1': 5.0, 'tau2': 16, 'k1': 0.00048, 'c1': 0.0024, 'tauplus': 22, \\\n",
    "                         'tauA': 33, 'k2': 560, 'dt': 1000/150.}\n",
    "        filename = './data/File078Moff.txt'\n",
    "        celltype = 'off-cell'\n",
    "        fitparameters_s = fitparameters.copy()\n",
    "        fitparameters_s['rsurr'] = fitparameters_s['rcenter']\n",
    "        retinal_reliability = 0.89**2 # From 78 MC off cell\n",
    "\n",
    "\n",
    "    # full model\n",
    "    rates_center, rates_surround = filter_video_MC('../stimuli/1x10_256.mpg', 9750, fitparameters)\n",
    "    RpcF = get_MC_response(rates_center, rates_surround, fitparameters, type = celltype)\n",
    "\n",
    "    # no-surround model\n",
    "    rates_center, rates_surround = filter_video_MC('../stimuli/1x10_256.mpg', 9750, fitparameters_s)\n",
    "    RpcS = get_MC_response(rates_center, rates_surround, fitparameters_s, type = celltype)\n",
    "    \n",
    "    make_plots(cell, RpcF, RpcS, filename, dt = fitparameters['dt'], shift = 2, column = 0, retinal_reliability = retinal_reliability)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Blue ON "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = \"178Bon\"\n",
    "\n",
    "fitparameters = {'xoff': -7.6, 'yoff': 9.5, 'rcenterB': 10.0, 'rcenterY': 14.0,\\\n",
    "                 'tau1': 10.2, 'tau2': 86, 'delay': 29, 'beta': 0.0, 'alpha': 4.9,  \\\n",
    "                 'tau3': 280, 'kL': 0.01, 'kM': 0.051, 'kB': 0.12, 'offset': 0.025, 'dt': 1000/150.}\n",
    "\n",
    "filename = './data/File178Bon.txt'\n",
    "fitparameters_s = fitparameters.copy()\n",
    "fitparameters_s['rcenterY'] = fitparameters_s['rcenterB']\n",
    "retinal_reliability = 0.91**2 # From barry 178 Bon cell\n",
    "\n",
    "\n",
    "# full model\n",
    "lum_signalL, lum_signalM, lum_signalB = filter_video_BO('../stimuli/1x10_256.mpg', 9750, fitparameters)\n",
    "RpcF = get_BO_response(lum_signalB, lum_signalM, lum_signalL, fitparameters)\n",
    "\n",
    "# no-surround model\n",
    "lum_signalL, lum_signalM, lum_signalB = filter_video_BO('../stimuli/1x10_256.mpg', 9750, fitparameters_s)\n",
    "RpcS = get_BO_response(lum_signalB, lum_signalM, lum_signalL, fitparameters_s)\n",
    "\n",
    "make_plots(cell, RpcF, RpcS, filename, dt = fitparameters['dt'], shift = 2, column = 0, retinal_reliability = retinal_reliability)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}