Neural_Activity_Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Activity Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relevant Packages and Utility Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Recompiling stale cache file C:\\Users\\Daniel Keitley\\.julia\\compiled\\v1.2\\ImageFiltering\\MN5uU.ji for ImageFiltering [6a3955dd-da59-5b1f-98d4-e7296123deb5]\n",
      "└ @ Base loading.jl:1240\n",
      "┌ Info: Precompiling PyPlot [d330b81b-6aea-500a-939a-2ce795aea3ee]\n",
      "└ @ Base loading.jl:1242\n",
      "ERROR: LoadError: LoadError: PyCall not properly installed. Please run Pkg.build(\"PyCall\")\n",
      "Stacktrace:\n",
      " [1] error(::String) at .\\error.jl:33\n",
      " [2] top-level scope at C:\\Users\\Daniel Keitley\\.julia\\packages\\PyCall\\eAqOC\\src\\startup.jl:44\n",
      " [3] include at .\\boot.jl:328 [inlined]\n",
      " [4] include_relative(::Module, ::String) at .\\loading.jl:1094\n",
      " [5] include at .\\Base.jl:31 [inlined]\n",
      " [6] include(::String) at C:\\Users\\Daniel Keitley\\.julia\\packages\\PyCall\\eAqOC\\src\\PyCall.jl:3\n",
      " [7] top-level scope at C:\\Users\\Daniel Keitley\\.julia\\packages\\PyCall\\eAqOC\\src\\PyCall.jl:37\n",
      " [8] include at .\\boot.jl:328 [inlined]\n",
      " [9] include_relative(::Module, ::String) at .\\loading.jl:1094\n",
      " [10] include(::Module, ::String) at .\\Base.jl:31\n",
      " [11] top-level scope at none:2\n",
      " [12] eval at .\\boot.jl:330 [inlined]\n",
      " [13] eval(::Expr) at .\\client.jl:432\n",
      " [14] top-level scope at .\\none:3\n",
      "in expression starting at C:\\Users\\Daniel Keitley\\.julia\\packages\\PyCall\\eAqOC\\src\\startup.jl:41\n",
      "in expression starting at C:\\Users\\Daniel Keitley\\.julia\\packages\\PyCall\\eAqOC\\src\\PyCall.jl:37\n",
      "ERROR: LoadError: Failed to precompile PyCall [438e738f-606a-5dbb-bf0a-cddfbfd45ab0] to C:\\Users\\Daniel Keitley\\.julia\\compiled\\v1.2\\PyCall\\GkzkC.ji.\n",
      "Stacktrace:\n",
      " [1] error(::String) at .\\error.jl:33\n",
      " [2] compilecache(::Base.PkgId, ::String) at .\\loading.jl:1253\n",
      " [3] _require(::Base.PkgId) at .\\loading.jl:1013\n",
      " [4] require(::Base.PkgId) at .\\loading.jl:911\n",
      " [5] require(::Module, ::Symbol) at .\\loading.jl:906\n",
      " [6] include at .\\boot.jl:328 [inlined]\n",
      " [7] include_relative(::Module, ::String) at .\\loading.jl:1094\n",
      " [8] include(::Module, ::String) at .\\Base.jl:31\n",
      " [9] top-level scope at none:2\n",
      " [10] eval at .\\boot.jl:330 [inlined]\n",
      " [11] eval(::Expr) at .\\client.jl:432\n",
      " [12] top-level scope at .\\none:3\n",
      "in expression starting at C:\\Users\\Daniel Keitley\\.julia\\packages\\PyPlot\\4wzW1\\src\\PyPlot.jl:4\n"
     ]
    },
    {
     "ename": "ErrorException",
     "evalue": "Failed to precompile PyPlot [d330b81b-6aea-500a-939a-2ce795aea3ee] to C:\\Users\\Daniel Keitley\\.julia\\compiled\\v1.2\\PyPlot\\oatAj.ji.",
     "output_type": "error",
     "traceback": [
      "Failed to precompile PyPlot [d330b81b-6aea-500a-939a-2ce795aea3ee] to C:\\Users\\Daniel Keitley\\.julia\\compiled\\v1.2\\PyPlot\\oatAj.ji.",
      "",
      "Stacktrace:",
      " [1] error(::String) at .\\error.jl:33",
      " [2] compilecache(::Base.PkgId, ::String) at .\\loading.jl:1253",
      " [3] _require(::Base.PkgId) at .\\loading.jl:1013",
      " [4] require(::Base.PkgId) at .\\loading.jl:911",
      " [5] require(::Module, ::Symbol) at .\\loading.jl:906",
      " [6] top-level scope at In[42]:1"
     ]
    }
   ],
   "source": [
    "using ImageFiltering, PyPlot, Plots, Plots.PlotMeasures;\n",
    "\n",
    "\n",
    "function report_progress(i::Int,iterations::Int,tstart::Float64)\n",
    "    \"\"\" Prints to screen the model progress. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    i: Int \n",
    "        Current iteration of the model.\n",
    "\n",
    "    iterations: Int \n",
    "        Number of iterations set to run. \n",
    "\n",
    "    tstart: Float64 \n",
    "        Time of starting the model. \n",
    "    \"\"\"\n",
    "    \n",
    "    progress = (i*100)/iterations\n",
    "    t = time()\n",
    "    tpassed = (t-tstart)/60\n",
    "    time2go = ((tpassed/progress)*100) - tpassed\n",
    "\n",
    "    print(\"$(@sprintf(\"%.3f\",progress))% completed.  Time passed: $(@sprintf(\"%.3f\",tpassed)) mins.    Expected time left: $(@sprintf(\"%.3f\",time2go))  \\r\")\n",
    "    flush(STDOUT)\n",
    "\n",
    "end\n",
    "\n",
    "\n",
    "function gen_pairs(ni,nj)\n",
    "    \"\"\" Generates all possible pairs of adjacent retinal neurons. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ni: Int\n",
    "        Number of retinal cells along nasal-temporal axis.\n",
    "    nj: Int\n",
    "        Number of retinal cells along dorsal-ventral axis.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    num_pairs: Int\n",
    "        Total number of possible pairs of adjacent neurons\n",
    "        in the presynaptic sheet.\n",
    "\n",
    "    pairs: 2D Array\n",
    "        Takes format pairs[pair_num,:] = [i1,j1,i2,j2] where i1,i2 and j1,j2 are\n",
    "        retinal coordinates in the nasal-temporal and dorsal-ventral axes\n",
    "        respectively. \n",
    "    \"\"\"\n",
    "    \n",
    "    num_pairs = 2*(ni-1)*(nj-1) + (ni-1) + (nj-1)\n",
    "    pairs = zeros(Int,(num_pairs,4))\n",
    "    count = 1 \n",
    "    for i=1:ni, j=1:nj\n",
    "        if i!=ni\n",
    "            pairs[count,:] = [i,j,i+1,j]\n",
    "            count += 1\n",
    "        end\n",
    "\n",
    "        if j!=nj\n",
    "            pairs[count,:] = [i,j,i,j+1]\n",
    "            count += 1\n",
    "        end\n",
    "    end\n",
    "\n",
    "    return (num_pairs,pairs)\n",
    "\n",
    "end\n",
    "\n",
    "\n",
    "function get_tec_pms(nk,nl)\n",
    "    \"\"\" Obtains coordinates of middle 4 tectal cells. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    nk: Int\n",
    "        Number of tectal cells along posterior-anterior axis.\n",
    "\n",
    "    nl: Int\n",
    "        Number of tectal cells along lateral-medial axis. \n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    pms: Array{Array{Int64,1},1}\n",
    "        k,l coordinates of the middle 4 tectal neurons. \n",
    "    \"\"\" \n",
    "    \n",
    "    kmid1 = fld(nk,2) \n",
    "    kmid2 = kmid1+1\n",
    "\n",
    "    lmid1 = fld(nl,2) \n",
    "    lmid2 = lmid1 + 1\n",
    "    \n",
    "    pms = [[kmid1,lmid1],[kmid1,lmid2],[kmid2,lmid1],[kmid2,lmid2]]\n",
    "\n",
    "    return pms\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: centered not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: centered not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[2]:18"
     ]
    }
   ],
   "source": [
    "const FRAME_RATE = 500\n",
    "const dt = 1\n",
    "\n",
    "const INIT_SYN_MEAN = 2.5                                 # S\n",
    "const INIT_SYN_SD = 0.14\n",
    "\n",
    "const POLARITY_MULT = 5\n",
    "const ACTIVITY_THRESH = 10                                # theta\n",
    "const MEM_TIME_CONST = 0.5                                # alpha\n",
    "const STATIONARITY_CONST = 0.005\n",
    "const SYN_MOD_THRESH = 2                                  # epsilon\n",
    "const HEBB_CONST = 0.0016                                 # h\n",
    "const INPUT_CLUSTER_SIZE = 2                              # c\n",
    "\n",
    "const beta = 0.05\n",
    "const gamma = 0.025\n",
    "const delta = -0.06\n",
    "const DIST_FILTER = centered([0 0 0 delta 0 0 0;\n",
    "            0 0 delta gamma delta 0 0;\n",
    "        0 delta gamma beta gamma delta 0;\n",
    "            delta gamma beta 0 beta gamma delta;\n",
    "            0 delta gamma beta gamma delta 0;\n",
    "                0 0 delta gamma delta 0 0;\n",
    "                    0 0 0 delta 0 0 0]);\n",
    "const offsets = [-3,-2,-1,0,1,2,3];\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialise Synaptic Weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "init_syns (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function init_syns(ni,nj,nk,nl,tec_pms)\n",
    "    \"\"\" Initialises synaptic weights according to a normal distribution.\n",
    "    Polarity markers are added to the middle 4 cells of each sheet.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ni: Int\n",
    "        Number of retinal neurons along nasal-temporal axis. \n",
    "\n",
    "    nj: Int\n",
    "        Number of retinal neurons along dorsal-ventral axis.\n",
    "\n",
    "    nk: Int\n",
    "        Number of tectal cells along posterior-anterior axis.\n",
    "\n",
    "    nl: Int\n",
    "        Number of tectal cells along lateral-medial axis.\n",
    "\n",
    "    tec_pms: Array{Array{Int64,1},1}\n",
    "        k,l coordinates of tectal polarity markers. \n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    syn_weights: Array{Float64,4}\n",
    "        Synaptic weights of connections between each retinal and tectal neuron.\n",
    "    \"\"\"\n",
    "    \n",
    "    syn_weights = randn(ni,nj,nk,nl)*INIT_SYN_SD + INIT_SYN_MEAN\n",
    "\n",
    "    imid1 =  fld(ni,2)\n",
    "    imid2 = imid1 + 1\n",
    "\n",
    "    jmid1 = fld(nj,2)\n",
    "    jmid2 = jmid1 + 1\n",
    "\n",
    "    tmid1 = tec_pms[1]\n",
    "    k1 = tmid1[1]\n",
    "    l1 = tmid1[2]\n",
    "\n",
    "    tmid2 = tec_pms[2]\n",
    "    k2 = tmid2[1]\n",
    "    l2 =tmid2[2]\n",
    "\n",
    "    tmid3 = tec_pms[3]\n",
    "    k3 = tmid3[1]\n",
    "    l3 = tmid3[2]\n",
    "\n",
    "    tmid4 = tec_pms[4]\n",
    "    k4 = tmid4[1]\n",
    "    l4 = tmid4[2]\n",
    "\n",
    "    syn_weights[imid1,jmid1,k1,l1] *= POLARITY_MULT\n",
    "    syn_weights[imid1,jmid2,k2,l2] *= POLARITY_MULT\n",
    "    syn_weights[imid2,jmid1,k3,l3] *= POLARITY_MULT\n",
    "    syn_weights[imid2,jmid2,k4,l4] *= POLARITY_MULT\n",
    "\n",
    "    syn_weights = norm_syns(syn_weights,ni,nj,nk,nl)\n",
    "\n",
    "    return syn_weights\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constructing model object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "syntax: extra token \"ActivityModel\" after end of expression",
     "output_type": "error",
     "traceback": [
      "syntax: extra token \"ActivityModel\" after end of expression",
      ""
     ]
    }
   ],
   "source": [
    "type ActivityModel\n",
    "    ret_dim1::Int\n",
    "    ret_dim2::Int\n",
    "    tec_dim1::Int\n",
    "    tec_dim2::Int\n",
    "    tec_pms::Array{Array{Int64,1},1}\n",
    "    pairs::Array{Int,2}\n",
    "    num_pairs::Int\n",
    "    syn_weights::Array{Float64,4}\n",
    "    rf_centres::Array{Float64,3}\n",
    "    iter_weights::Array{Float64,5}\n",
    "    iter_rfs::Array{Float64,4}\n",
    "    cur_pot::Array{Float64,2}\n",
    "    activity::Array{Float64,2}\n",
    "    conv::Array{Float64,2}\n",
    "    pot_change::Array{Float64,2}\n",
    "    iinds::Array{Int,1}\n",
    "    jinds::Array{Int,1}\n",
    "\n",
    "    ## Constructor\n",
    "    function ActivityModel(ni,nj,nk,nl,tec_pms)\n",
    "        s = init_syns(ni, nj, nk, nl,tec_pms)\n",
    "        (num_pairs,pairs) = gen_pairs(ni,nj)\n",
    "        rf = zeros(nk, nl,2)\n",
    "        iw = zeros(1,ni,nj,nk,nl)\n",
    "        irf = zeros(1,nk,nl,2)\n",
    "        cur_pot = zeros(nk,nl)\n",
    "        activity = zeros(nk,nl)\n",
    "        conv = zeros(nk,nl)\n",
    "        pot_change = zeros(nk,nl)\n",
    "        iinds = repeat([ i for i in 1:ni],inner=nj)\n",
    "        jinds = repeat([ j for j in 1:nj],outer=ni)\n",
    "\n",
    "        new(ni, nj, nk, nl, tec_pms,pairs,num_pairs, s, rf,iw,irf,\n",
    "            cur_pot,activity,conv,pot_change,iinds,jinds)\n",
    "    end\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Activate adjacent retinal cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "activate_ret (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function activate_ret(pairs,num_pairs)\n",
    "    \"\"\" Chooses a random pair of adjacent retinal neurons to activate.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    pairs: Array{Int,2}\n",
    "        i,j coordinates of all adjacent pairs of retinal neurons. \n",
    "        Takes the form pairs[pair_num,:] = [i1,j1,i2,j2].\n",
    "    num_pairs: int \n",
    "        Number of distinct adjacent pairs of retinal neurons.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    pair: Array{Int,1}\n",
    "        i,j coordinates of the chosen pair of retinal neurons. \n",
    "    \"\"\"\n",
    "    \n",
    "    pair = Array{Int}(1)\n",
    "    pair  = pairs[rand(1:num_pairs),:]\n",
    "    return pair\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Depolarise tectal cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[6]:1"
     ]
    }
   ],
   "source": [
    "function activate_tec(model::ActivityModel,pair)\n",
    "    \"\"\" Computes linear sum of activated retinal neurons and their\n",
    "    synaptic weights. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component.\n",
    "\n",
    "    pair: Array{Int,1} \n",
    "        i,j, coordinates of activated retinal neurons. \n",
    "    \"\"\"\n",
    "    \n",
    "    i1 = pair[1]\n",
    "    j1 = pair[2]\n",
    "    i2 = pair[3]\n",
    "    j2 = pair[4]\n",
    "\n",
    "    model.cur_pot = model.syn_weights[i1,j1,:,:] + model.syn_weights[i2,j2,:,:]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Threshold tectal depolarisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "thresh (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function thresh(cur_pot)\n",
    "    \"\"\" Applies threshold to tectal membrane potentials. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    cur_pot: Array{Float64,2}\n",
    "        Membrane potential in each tectal neuron. \n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    cur_pot: Array{Float64,2}\n",
    "        Post-threshold activity of each tectal neuron. \n",
    "    \"\"\"\n",
    "    \n",
    "    cur_pot[find(x -> x < ACTIVITY_THRESH, cur_pot)] = 0\n",
    "    cur_pot[find(x -> x > ACTIVITY_THRESH, cur_pot)] -= ACTIVITY_THRESH\n",
    "    return cur_pot\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Apply tectal excitation/inhibition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[8]:1"
     ]
    }
   ],
   "source": [
    "function tec_feedback(model::ActivityModel)\n",
    "    \"\"\" Computes excitation and inhibition effects from neighbouring tectal\n",
    "    neurons. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component.  \n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    feedback: Array{Float64,2} \n",
    "        Convolution of tectal activity with a Manhatten distance filter. \n",
    "        Represents contributions of neighbouring tectal neurons to \n",
    "        membrane potential. \n",
    "    \"\"\"\n",
    "    \n",
    "    nk = model.tec_dim1\n",
    "    nl = model.tec_dim2\n",
    "    feedback = zeros(Float64,nk,nl)\n",
    "    \n",
    "    for k in 1:nk, l in 1:nl, koff in offsets, loff in offsets\n",
    "        kind = k+koff\n",
    "        lind = l+loff\n",
    "        if kind>0 && lind >0 && kind <= nk && lind <= nl\n",
    "            feedback[k,l] += model.activity[kind,lind]*DIST_FILTER[koff,loff]\n",
    "        end        \n",
    "    end\n",
    "\n",
    "    return feedback\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Update membrane potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[9]:1"
     ]
    }
   ],
   "source": [
    "function update_potential(model::ActivityModel)\n",
    "    \"\"\" Simulates differential equations, using Euler's method, to update\n",
    "    tectal membrane potentials (see Willshaw, von der Malsburg, 1976). \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component.\n",
    "    \"\"\"\n",
    "    \n",
    "    init_pot = Array{Float64}(2)\n",
    "    init_pot = copy(model.cur_pot)\n",
    "    mean_change = 1\n",
    "    new_mean = mean(model.cur_pot)\n",
    "    old_mean = 0\n",
    "    while mean_change > (STATIONARITY_CONST*old_mean)\n",
    "        old_mean = new_mean\n",
    "        model.activity = thresh(copy(model.cur_pot))\n",
    "        model.conv = tec_feedback(model)\n",
    "        model.pot_change = (init_pot + model.conv) - (MEM_TIME_CONST*model.cur_pot) \n",
    "        model.cur_pot += model.pot_change*dt\n",
    "        new_mean = mean(model.cur_pot)\n",
    "        mean_change = abs(new_mean-old_mean)\n",
    "    end\n",
    "\n",
    "    model.activity = thresh(model.cur_pot)\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Apply Hebbian learning rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[10]:1"
     ]
    }
   ],
   "source": [
    "function update_syns(model::ActivityModel,pair)\n",
    "    \"\"\" Applies Hebbian learning rule to connections between activated \n",
    "    retinal neurons and strongly excitated tectal neurons. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component.\n",
    "\n",
    "    pair: Array{Int,1} \n",
    "        i,j, coordinates of activated retinal neurons.  \n",
    "    \"\"\"\n",
    "    \n",
    "    i1 = pair[1]\n",
    "    j1 = pair[2]\n",
    "    i2 = pair[3]\n",
    "    j2 = pair[4]\n",
    "\n",
    "    linds = find(x -> x > SYN_MOD_THRESH,model.activity)\n",
    "\n",
    "    if length(linds)!=0\n",
    "        for lind in linds\n",
    "            (k,l) = ind2sub(model.activity,lind)\n",
    "            model.syn_weights[i1,j1,k,l] += HEBB_CONST*model.activity[k,l]\n",
    "            model.syn_weights[i2,j2,k,l] += HEBB_CONST*model.activity[k,l]\n",
    "        end\n",
    "    end\n",
    "\n",
    "    model.syn_weights = norm_syns(model.syn_weights,model.ret_dim1,\n",
    "        model.ret_dim2,model.tec_dim1,model.tec_dim2)\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Normalise synaptic weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "norm_syns (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function norm_syns(syn_weights,ni,nj,nk,nl)\n",
    "    \"\"\" Normalises the synaptic weights so that the mean weight \n",
    "    contacting each tectal neuron remains constant. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    syn_weights: Array{Float64,4} \n",
    "        Synaptic weights of each connection between retinal and tectal neurons.\n",
    "\n",
    "    ni: Int\n",
    "        Number of retinal neurons along nasal-temporal axis. \n",
    "\n",
    "    nj: Int\n",
    "        Number of retinal neurons along dorsal-ventral axis.\n",
    "\n",
    "    nk: Int\n",
    "        Number of tectal cells along posterior-anterior axis.\n",
    "\n",
    "    nl: Int\n",
    "        Number of tectal cells along lateral-medial axis.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    syn_weights: Array{Float64,4}\n",
    "        Normalised synaptic weights. \n",
    "    \"\"\"\n",
    "    \n",
    "    reshaped_weights = reshape(syn_weights,(ni*nj,nk*nl))\n",
    "    mu = mean(reshaped_weights,1)\n",
    "    norm_weights = (reshaped_weights*INIT_SYN_MEAN)./mu\n",
    "    syn_weights = reshape(norm_weights,(ni,nj,nk,nl))\n",
    "    return syn_weights\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Putting it all together..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[12]:1"
     ]
    }
   ],
   "source": [
    "function iteration(model::ActivityModel)\n",
    "    \"\"\" Simulates one cycle of the model simulation. \n",
    "\n",
    "    A simulation cycle consists of:\n",
    "    1) Exciting a pair of retinal neurons at random.\n",
    "    2) Computing the subseqent activation in the tectum.\n",
    "    3) Simulating tectal interactions.\n",
    "    4) Updating synaptic weights. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component.\n",
    "    \"\"\"\n",
    "    \n",
    "    pair = activate_ret(model.pairs,model.num_pairs)\n",
    "    activate_tec(model,pair)\n",
    "    update_potential(model)\n",
    "    update_syns(model,pair)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determining tectal receptive fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[13]:1"
     ]
    }
   ],
   "source": [
    "function calc_rf_centres(model::ActivityModel)\n",
    "    \"\"\" Computes retinal coordinates for the receptive field centre of each\n",
    "    tectal neuron. For each tectal neuron, the coordinates of the receptive\n",
    "    field centre is calculated as a sum over all retinal positions \n",
    "    weighted by the corresponding synatpic weights. The weighted\n",
    "    sum is then normalised by dividing over the sum of synaptic weights \n",
    "    arriving at the tectal neuron being considered.  \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel \n",
    "        ActivityModel object containing current values of each model \n",
    "        component.\n",
    "    \"\"\"\n",
    "    \n",
    "    ni = model.ret_dim1\n",
    "    nj = model.ret_dim2\n",
    "    nk = model.tec_dim1\n",
    "    nl = model.tec_dim2\n",
    "\n",
    "    rshpe_syns = reshape(model.syn_weights,(ni*nj,nk*nl))\n",
    "    rfdenom = sum(rshpe_syns,1)\n",
    "\n",
    "    rfnume = sum((rshpe_syns.*model.iinds),1)\n",
    "    irf = reshape(rfnume./rfdenom,(nk,nl))\n",
    "\n",
    "    rfnume = sum((rshpe_syns.*model.jinds),1)\n",
    "    jrf = reshape(rfnume./rfdenom,(nk,nl))\n",
    "\n",
    "    model.rf_centres[:,:,1] = irf\n",
    "    model.rf_centres[:,:,2] = jrf\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running Model Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[14]:1"
     ]
    }
   ],
   "source": [
    "function run(model::ActivityModel,iterations::Int)\n",
    "    \"\"\" Runs simulation of Activity Model.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component.\n",
    "        \n",
    "    iterations: Int\n",
    "        Number of iterations to iterate simulation cycle. \n",
    "    \"\"\"\n",
    "    \n",
    "    ni = model.ret_dim1\n",
    "    nj = model.ret_dim2\n",
    "    nk = model.tec_dim1\n",
    "    nl = model.tec_dim2\n",
    "\n",
    "    model.iter_weights = zeros(floor(Int64,iterations/FRAME_RATE)+1,\n",
    "        ni,nj,nk,nl)\n",
    "    model.iter_rfs = zeros(floor(Int64,iterations/FRAME_RATE)+1,\n",
    "        nk,nl,2)\n",
    "\n",
    "    tstart = time()\n",
    "\n",
    "    for i in 1:iterations\n",
    "        iteration(model)\n",
    "\n",
    "        if i % FRAME_RATE == 0\n",
    "            index = floor(Int64,i/FRAME_RATE)\n",
    "            model.iter_weights[index,:,:,:,:] = model.syn_weights\n",
    "            calc_rf_centres(model)\n",
    "            model.iter_rfs[index,:,:,:] = model.rf_centres\n",
    "            report_progress(i,iterations,tstart)\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualising Retinotopic Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: ActivityModel not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: ActivityModel not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[15]:1"
     ]
    }
   ],
   "source": [
    "function plot_iteration(model::ActivityModel,iteration::Int,title::String)\n",
    "    \"\"\" Plots a lattice plot representing the receptive fields of \n",
    "    each tectal neuron at a given iteration of the model simulation. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component. \n",
    "\n",
    "    iteration: int \n",
    "        Iteration of model simulation to display receptive fields. \n",
    "    \n",
    "    title: String\n",
    "        Text to display as subplot title. \n",
    "    \"\"\"\n",
    "    \n",
    "    rec_fields = model.iter_rfs[iteration,:,:,:]\n",
    "    ax = PyPlot.gca()\n",
    "    ax[:set_xlim](1,model.ret_dim1)\n",
    "    ax[:set_ylim](1,model.ret_dim2)\n",
    "    ax[:set_aspect](\"equal\")\n",
    "    PyPlot.grid(\"on\")\n",
    "    \n",
    "    for k in 1:model.tec_dim1\n",
    "        PyPlot.plot(rec_fields[k,:,1],rec_fields[k,:,2],c=\"k\")\n",
    "    end\n",
    "\n",
    "\n",
    "    for l in 1:model.tec_dim2\n",
    "        PyPlot.plot(rec_fields[:,l,1],rec_fields[:,l,2],c=\"k\")\n",
    "    end\n",
    "    \n",
    "    PyPlot.title(title)\n",
    "    PyPlot.xlabel(\"Nasal-Temporal Position\")\n",
    "    PyPlot.ylabel(\"Dorsal-Ventral Position\")\n",
    "  \n",
    "end\n",
    "\n",
    "\n",
    "\n",
    "function plot_map(model::ActivityModel)\n",
    "    \"\"\" Sets up lattice plot figure and displays the receptive field \n",
    "    centres at 4 equally spaced iterations of the model simulation. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component. \n",
    "    \"\"\"\n",
    "    fig = figure()\n",
    "    iterations = (length(model.iter_rfs[:,1,1,1])-1)\n",
    "    PyPlot.subplot(2,2,1)\n",
    "    plot_iteration(model,1,\"Iteration: 1\")\n",
    "    \n",
    "    PyPlot.subplot(2,2,2)\n",
    "    it = fld(iterations,3)\n",
    "    plot_iteration(model,it,\"Iteration: $(it*FRAME_RATE)\")\n",
    "    \n",
    "    PyPlot.subplot(2,2,3)\n",
    "    it = fld(2*iterations,3)\n",
    "    plot_iteration(model,it,\"Iteration: $(it*FRAME_RATE)\")\n",
    "    \n",
    "    PyPlot.subplot(2,2,4)\n",
    "    it = iterations-1\n",
    "    plot_iteration(model,it,\"Iteration: $(it*FRAME_RATE)\")\n",
    "    \n",
    "    PyPlot.tight_layout()\n",
    "\n",
    "end\n",
    "\n",
    "\n",
    "function animate_map(model::ActivityModel, frames::Int,fps::Int)\n",
    "    \"\"\" Displays an mp4 animated movie of the model simulations.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model: ActivityModel\n",
    "        ActivityModel object containing current values of each model \n",
    "        component. \n",
    "    \n",
    "    frames: Int\n",
    "        Number of frames to animate. \n",
    "    \n",
    "    fps: Int\n",
    "        Frame rate of video produced. \n",
    "    \"\"\"\n",
    "    \n",
    "    anim = Animation()\n",
    "\n",
    "    for i in 1:frames\n",
    "        rec_fields = model.iter_rfs[i,:,:,:]\n",
    "        Plots.plot()\n",
    "\n",
    "        for k in 1:model.tec_dim1\n",
    "            Plots.plot!(rec_fields[k,:,1],rec_fields[k,:,2],legend=false,c=\"black\",\n",
    "                xlims = (1,model.ret_dim1),ylims=(1,model.ret_dim2),aspect_ratio=1,title=\"Iteration: $(i)\",\n",
    "                xlabel=\"Nasal-Temporal Position\",ylabel=\"Dorsal-Ventral Position\")\n",
    "        end\n",
    "\n",
    "        for l in 1:model.tec_dim2\n",
    "            Plots.plot!(rec_fields[:,l,1],rec_fields[:,l,2],legend=false,c=\"black\",\n",
    "                xlims = (1,model.ret_dim1),ylims=(1,model.ret_dim2),aspect_ratio=1,title=\"Iteration: $(i)\",\n",
    "            xlabel=\"Nasal-Temporal Position\",ylabel=\"Dorsal-Ventral Position\")\n",
    "        end\n",
    "\n",
    "        frame(anim)\n",
    "    end\n",
    "\n",
    "    mp4(anim,fps=fps)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: get_tec_pms not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: get_tec_pms not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[16]:1"
     ]
    }
   ],
   "source": [
    "am1 = ActivityModel(6,6,6,6,get_tec_pms(6,6));\n",
    "run(am1,60000);\n",
    "plot_map(am1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Systems Matching"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6x6 | 9x6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: get_tec_pms not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: get_tec_pms not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[17]:1"
     ]
    }
   ],
   "source": [
    "am2 = ActivityModel(6,6,9,6,get_tec_pms(9,6));\n",
    "run(am2,100000);\n",
    "plot_map(am2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8x8 | 6x6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: get_tec_pms not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: get_tec_pms not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[18]:1"
     ]
    }
   ],
   "source": [
    "am3 = ActivityModel(8,8,6,6,get_tec_pms(6,6));\n",
    "run(am3,80000);\n",
    "plot_map(am3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8x8 | 9x6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: get_tec_pms not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: get_tec_pms not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[19]:1"
     ]
    }
   ],
   "source": [
    "am4 = ActivityModel(8,8,9,6,get_tec_pms(9,6));\n",
    "run(am4,120000);\n",
    "plot_map(am4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing Polarity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "change_polarity (generic function with 1 method)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function change_polarity(perm,nk,nl)\n",
    "    tec_pms = get_tec_pms(nk,nl)\n",
    "    cop = copy(tec_pms)\n",
    "    display()\n",
    "    for n in 1:4\n",
    "        tec_pms[n] = cop[perm[n]] \n",
    "    end\n",
    "    \n",
    "    return tec_pms\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: get_tec_pms not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: get_tec_pms not defined",
      "",
      "Stacktrace:",
      " [1] change_polarity(::Array{Int64,1}, ::Int64, ::Int64) at .\\In[20]:2",
      " [2] top-level scope at In[21]:1"
     ]
    }
   ],
   "source": [
    "tec_pms1 = change_polarity([2,4,1,3],9,6)\n",
    "am5 = ActivityModel(8,8,9,6,tec_pms1);\n",
    "run(am5,120000);\n",
    "plot_map(am5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Twist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: get_tec_pms not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: get_tec_pms not defined",
      "",
      "Stacktrace:",
      " [1] change_polarity(::Array{Int64,1}, ::Int64, ::Int64) at .\\In[20]:2",
      " [2] top-level scope at In[22]:1"
     ]
    }
   ],
   "source": [
    "tec_pms2 = change_polarity([1,4,3,2],9,6)\n",
    "am6 = ActivityModel(8,8,9,6,tec_pms2);\n",
    "run(am6,120000);\n",
    "plot_map(am6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}