#!/usr/bin/env python
# coding: utf-8

# ## SYDE556/750 Assignment 1: Representation in Populations of Neurons
# 
# - Due Date: January 28th at midnight
# - Total marks: 20 (20% of final grade)
# - Late penalty: 1 mark per day
# - It is recommended that you use Python.
# - *Do not use or refer to any code from Nengo*

# ## 1) Representation of Scalars

# ### 1.1) Basic encoding and decoding
# 
# Write a program that implements a neural representation of a scalar value $x$.  For the neuron model, use a rectified linear neuron model ($a=max(J,0)$).  Choose the maximum firing rates randomly (uniformly distributed between 100Hz and 200Hz at x=1), and choose the x-intercepts randomly (uniformly distributed between -0.95 and 0.95).  Use those values to compute the corresponding $\alpha$ and $J^{bias}$ parameters for each neuron.  The encoders $e$ are randomly chosen and are either +1 or -1 for each neuron.  Go through the following steps:
# 
# 
# 

# <ol style="list-style:lower-alpha outside none;" start="1">
# <li>[1 mark] Plot the neuron responses $a_i$ for 16 randomly generated neurons. (See Figure 2.4 in the book for an example, but with a different neuron model and a different range of maximum firing rates).  
# <ul><li>Since you can't compute this for every possible $x$ value between -1 and 1, sample the x-axis with $dx=0.05$.  Use this sampling throughout this question)</li></ul>
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="2">
# <li> [1 mark] Compute the optimal decoders $d_i$ for those 16 neurons (as shown in class).  Report their values.
# <ul><li>The easiest way to compute $d$ is to use the matrix notation mentioned in the course notes.  $A$ is the matrix of neuron activities (the same thing used to generate the plot in 1.1a).
# </li></ul>
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="3">
# <li> [1 mark] Compute and plot $\hat{x}=\sum_i d_i a_i$.  Overlay on the plot the line $y=x$.  (See Figure 2.7 for an example).  Make a separate plot of $x-\hat{x}$ to see what the error looks like.  Report the Root Mean Squared Error value.
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="4">
# <li> [1 mark] Now try decoding under noise.  Add random normally distributed noise to $a$ and decode again. The noise is a random variable with mean 0 and standard deviation of 0.2 times the maximum firing rate of all the neurons.  Resample this variable for every different $x$ value for every different neuron.  Create all the same plots as in part c).  Report the Root Mean Squared Error value.
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="5">
# <li>[1 mark] Recompute the decoders $d_i$ taking noise into account (as shown in class).  Show how these decoders behave when decoding both with and without noise added to $a$ by making the same plots as in c) and d).  Report the RMSE for all cases.
# <ul><li>As in the previous question, $\sigma$ is 0.2 times the maximum firing rate of all the neurons.
# </li></ul>
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="6">
# <li> [1 mark] Show a 2x2 table of the four RMSE values reported in parts c), d), and e).  This should show the effects of adding noise and whether or not the decoders $d$ are computed taking noise into account.  Write a few sentences commenting on what the table shows.
# </li>
# </ol>

# ### 1.2) Exploring sources of error
# 
# Use the program you wrote in 1.1 to examine the sources of error in the representation.

# <ol style="list-style:lower-alpha outside none;" start="1">
# <li>[2 marks] Plot the error due to distortion $E_{dist}$ and the error due to noise $E_{noise}$ as a function of $N$, the number of neurons. Use the equation with those two parts as your method (2.9 in the book). Generate two different loglog plots (one for each type of error) with $N$ values of [4, 8, 16, 32, 64, 128, 256, 512] (and more, if you would like).  For each $N$ value, do at least 5 runs and average the results.  For each run, different $\alpha$, $J^{bias}$, and $e$ values should be generated for each neuron.  Compute $d$ under noise, with $\sigma$ equal to 0.1 times the maximum firing rate.  Show visually that the errors are proportional to $1/N$ or $1/N^2$ (see figure 2.6 in the book).
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="2">
# <li>[1 mark] Repeat part a) with $\sigma$ equal to 0.01 times the maximum firing rate.
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="3">
# <li> [1 mark] What does the difference between the graphs in a) and b) tell us about the sources of error in neural populations?
# </li>
# </ol>

# ### 1.3) Leaky Integrate-and-Fire neurons
# 
# Change the code to use the LIF neuron model: 
# 
# $$
# a_i = \begin{cases}
#     {1 \over {\tau_{ref}-\tau_{RC}ln(1-{1 \over J})}} &\mbox{if } J>1 \\ 
#     0 &\mbox{otherwise} 
#     \end{cases}
# $$
# 

# <ol style="list-style:lower-alpha outside none;" start="1">
# <li> [1 mark] Generate the same plot as 1.1a).  Use $\tau_{ref}=0.002$s and $\tau_{RC}=0.02$s.
# <ul><li>Note that you will need to compute new $\alpha$ and $J^{bias}$ values that will achieve the desired tuning curves (uniform distribution of x-intercepts between -1 and 1, and maximum firing rates between 100Hz and 200Hz).  Since you know two points on the tuning curve (the x-intercept and the point where it hits maximum firing), this gives you 2 equations and 2 unknowns, so you can find $\alpha$ and $J^{bias}$ by substituting and rearranging.
# </li></ul>
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="2">
# <li> [2 marks] Generate the same four plots as 1.1e), and report the RMSE both with and without noise.
# </li>
# </ol>

# ## 2) Representation of Vectors

# ### 2.1) Vector tuning curves

# <ol style="list-style:lower-alpha outside none;" start="1">
# <li> [1 mark] Plot the tuning curve of an LIF neuron whose 2D preferred direction vector is at an angle of $\theta=-\pi/4$, has an x-intercept at the origin (0,0), and has a maximum firing rate of 100Hz.
# 
# <ul><li> Remember that $J=\alpha e \cdot x + J^{bias}$, and both $x$ and $e$ are 2D vectors.</li>
# <li> This is a 3D plot similar to figure 2.8a in the book.</li>
# <li> In the scalar case (that you did in question 1.1a), the maximum firing rate occurred when $x=1$ for neurons with $e=1$ and at $x=-1$ for neurons with $e=-1$.  Of course, if the graph in 1.1a was extended to $x>1$ (or $x<-1$), neurons would start firing faster than their maximum firing rate.  Similarly, here the "maximum firing rate" means the firing rate when $x=e$.  This should allow you to reuse your code from 1.3a) to compute $\alpha$ and $J^{bias}$ for a desired maximum firing rate and x-intercept.</li>
# <li> To generate 3D plots in MATLAB, see [here](http://www.mathworks.com/help/matlab/learn_matlab/creating-mesh-and-surface-plots.html)</li>
# <li> To generate 3D plots in Python, see [here](http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html)</li>
# </ul>
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="2">
# <li> [1 mark] Plot the tuning curve for the same neuron as in a), but only considering the points around the unit circle.  This will be similar to Figure 2.8b in the book.  Fit a curve of the form $Acos(B\theta+C)+D$ to the tuning curve and plot it as well.  What makes a cosine a good choice for this?  Why does it differ from the ideal curve?
# <ul><li> To do curve fitting in MATLAB, see [here](http://www.mathworks.com/help/optim/ug/lsqcurvefit.html).</li>
# <li> To do curve fitting in Python, see [here](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html).</li>
# </ul>
# </li>
# </ol>

# ### 2.2 Vector representation

# <ol style="list-style:lower-alpha outside none;" start="1">
# <li> [1 mark] Generate a set of 100 random unit vectors uniformly distributed around the unit circle.  These will be the encoders $e$ for 100 neurons.  Plot these vectors with a quiver or line plot (i.e. not just points, but lines/arrows to the points).
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="2">
# <li> [1 mark] Compute the optimal decoders.  Use LIF neurons with the same properties as in question 1.3.  When computing the decoders, take into account noise with $\sigma$ as 0.2 times the maximum firing rate.  Plot the decoders.  How do these decoding vectors compare to the encoding vectors?
# <ul><li> Note that the decoders will also be 2D vectors.</li>
# <li>In the scalar case, you used $x$ values between -1 and 1, with $dx=0.05$.  In this case, you can regularly tile the 2D $x$ values ([1, 1], [1, 0.95], ... [-1, -0.95], [-1, 1]).  Alternatively, you can just randomly choose 1600 different $x$ values to sample. </li>
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="3">
# <li> [1 mark] Generate 20 random $x$ values throughout the unit circle (i.e. with different directions and radiuses).  For each $x$ value, determine the neural activity $a$ for each of the 100 neurons.  Now decode these values (i.e. compute $\hat{x}$) using the decoders from part b).  Plot the original and decoded values on the same graph in different colours, and compute the RMSE.
# </li>
# </ol>

# <ol style="list-style:lower-alpha outside none;" start="4">
# <li> [2 marks] Repeat part c) but use the *encoders* as decoders.  This is what Georgopoulos used in his original approach to decoding information from populations of neurons.  Plot the decoded values this way and compute the RMSE.  In addition, recompute the RMSE in both cases you've done, but ignoring the magnitude of the decoded vector.  What are the relative merits of these two approaches to decoding?
# <ul><li>To ignore the magnitude of the vectors, normalize the length of the decoded vectors before computing the RMSE.</li>
# </ul>
# </li>
# </ol>

# In[ ]: