text:
# Lab 1: Independent Component Analysis

### Machine Learning 2 (2017/2018)

* The lab exercises should be made in groups of two people.
* The deadline is Thursday, April 19, 23:59.
* Assignment should be submitted through BlackBoard! Make sure to include your and your teammates' names with the submission.
* Attach the .IPYNB (IPython Notebook) file containing your code and answers. Naming of the file should be "studentid1\_studentid2\_lab#", for example, the attached file should be "12345\_12346\_lab1.ipynb". Only use underscores ("\_") to connect ids, otherwise the files cannot be parsed.

Notes on implementation:

* You should write your code and answers in an IPython Notebook: http://ipython.org/notebook.html. If you have problems, please ask.
* Use __one cell__ for code and markdown answers only!
    * Put all code in the cell with the ```# YOUR CODE HERE``` comment and overwrite the ```raise NotImplementedError()``` line.
    * For theoretical questions, put your solution using LaTeX style formatting in the YOUR ANSWER HERE cell.
* Among the first lines of your notebook should be "%pylab inline". This imports all required modules, and your plots will appear inline.
* Large parts of you notebook will be graded automatically. Therefore it is important that your notebook can be run completely without errors and within a reasonable time limit. To test your notebook before submission, select Kernel -> Restart \& Run All.

### Literature
In this assignment, we will implement the Independent Component Analysis algorithm as described in chapter 34 of David MacKay's book "Information Theory, Inference, and Learning Algorithms", which is freely available here:
http://www.inference.phy.cam.ac.uk/mackay/itila/book.html

Read the ICA chapter carefuly before you continue!

### Notation

$\mathbf{X}$ is the $M \times T$ data matrix, containing $M$ measurements at $T$ time steps.

$\mathbf{S}$ is the $S \times T$ source matrix, containing $S$ source signal values at $T$ time steps. We will assume $S = M$.

$\mathbf{A}$ is the mixing matrix. We have $\mathbf{X} = \mathbf{A S}$.

$\mathbf{W}$ is the matrix we aim to learn. It is the inverse of $\mathbf{A}$, up to indeterminacies (scaling and permutation of sources).

$\phi$ is an elementwise non-linearity or activation function, typically applied to elements of $\mathbf{W X}$.

### Code
In the following assignments, you can make use of the signal generators listed below.

```python
% pylab inline
import sys

assert sys.version_info[:3] >= (
3, 6, 0), "Make sure you have Python 3.6 installed"

# Signal generators
def sawtooth(x, period=0.2, amp=1.0, phase=0.):
    return (((x / period - phase - 0.5) % 1) - 0.5) * 2 * amp

def sine_wave(x, period=0.2, amp=1.0, phase=0.):
    return np.sin((x / period - phase) * 2 * np.pi) * amp

def square_wave(x, period=0.2, amp=1.0, phase=0.):
    return ((np.floor(2 * x / period - 2 * phase - 1) % 2 == 0).astype(
        float) - 0.5) * 2 * amp

def triangle_wave(x, period=0.2, amp=1.0, phase=0.):
    return (sawtooth(x, period, 1., phase) * square_wave(x, period, 1.,
                                                         phase) + 0.5) * 2 * amp

def random_nonsingular_matrix(d=2):
    """
    Generates a random nonsingular (invertible) matrix of shape d*d
    """
    epsilon = 0.1
    A = np.random.rand(d, d) #####change this function to randn to test with gaussian distribution
    while abs(np.linalg.det(A)) < epsilon:
        A = np.random.rand(d, d)
    return A

def plot_signals(X, title="Signals"):
    """
    Plot the signals contained in the rows of X.
    """
    figure()

    for i in range(X.shape[0]):
        ax = plt.subplot(X.shape[0], 1, i + 1)
        plot(X[i, :])
        ax.set_xticks([])
        ax.set_yticks([])
    plt.suptitle(title)

```

    Populating the interactive namespace from numpy and matplotlib

The following code generates some toy data to work with.

```python
# Generate data
num_sources = 5
signal_length = 500
t = linspace(0, 1, signal_length)
S = np.c_[sawtooth(t), sine_wave(t, 0.3), square_wave(t, 0.4), triangle_wave(t, 0.25), np.random.randn(t.size)].T
plot_signals(S)
```

### 1.1 Make mixtures (5 points)
Write a function `make_mixtures(S, A)' that takes a matrix of source signals $\mathbf{S}$ and a mixing matrix $\mathbf{A}$, and generates mixed signals $\mathbf{X}$.

```python
### 1.1 Make mixtures
def make_mixtures(S, A):
    # YOUR CODE HERE
    return np.dot(A, S)

```

```python
### Test your function
np.random.seed(42)
A = random_nonsingular_matrix(d=S.shape[0])
X = make_mixtures(S, A)
plot_signals(X, "Mixtures")

assert X.shape == (num_sources, signal_length), "The shape of your mixed signals is incorrect"

```

### 1.2 Histogram (5 points)
Write a function `plot_histograms(X)` that takes a data-matrix $\mathbf{X}$ and plots one histogram for each signal (row) in $\mathbf{X}$. You can use the `np.histogram()` (followed by `plot`) or `plt.hist()` function.

Plot histograms of the sources and the measurements.

```python
### 1.2 Histogram

import matplotlib.mlab as mlab
import math

def plot_histograms(X, title="Histogram"):
    mu = 0
    variance = 1
    sigma = math.sqrt(variance)
    x_ = np.linspace(mu - 2.5 * sigma, mu + 2.5 * sigma, 100)

    for index, x in enumerate(X):
        plt.subplot(1, len(X), index + 1)
        plt.hist(x)#, normed=True) # in order to compare with Gaussian distribution
        #plt.plot(x_, mlab.normpdf(x_, mu, sigma) , linewidth = 3) # in order to compare with Gaussian distribution

    plt.suptitle(title)
```

```python
plot_histograms(S, "Histograms of the Signals")
```

```python
plot_histograms(X, "Histograms of the Measurements")
```

ADD THEORETICAL ANSWER

Which of these distributions (sources or measurements) tend to look more like Gaussians? Can you think of an explanation for this phenomenon? Why is this important for ICA?

YOUR ANSWER HERE
The distributions of the measurements tend to look more like Gaussians. This is caused by (uniformly) randomly mixing the signals.
For the ICA it is important that the signals are non-Gaussian distributed in order to be able to separate them.

### 1.3 Implicit priors (20 points)
As explained in MacKay's book, an activation function $\phi$ used in the ICA learning algorithm corresponds to a prior distribution over sources. Specifically, $\phi(a) = \frac{d}{da} \ln p(a)$. For each of the following activation functions, *derive* the source distribution they correspond to.
$$\phi_0(a) = -\tanh(a)$$
$$\phi_1(a) = -a + \tanh(a)$$
$$\phi_2(a) = -a^3$$
$$\phi_3(a) = -\frac{6a}{a^2 + 5}$$

Give your answer without the normalizing constant, so an answer of the form $p(a) \propto \verb+[answer]+$ is ok.

YOUR ANSWER HERE

$$\phi(a) = \frac{d}{da} \ln p(a)$$
therefore:
$$p(a) \propto exp(\int \phi(a) da)$$

This gives us:
$$ p_0(a) \propto \frac{1}{cosh(a)} $$
$$ p_1(a) \propto exp(\frac{-a^2}{2})~ cosh(a) $$
$$ p_2(a) \propto exp(\frac{-a^4}{4})$$
$$ p_3(a) \propto \frac{1}{(a^2 + 5)^3} $$

```python
def phi_0(a):
    return -np.tanh(a)

def p_0(a):
    # YOUR CODE HERE
    return 1 / np.cosh(a)

```

```python
def phi_1(a):
    return -a + np.tanh(a)

def p_1(a):
    # YOUR CODE HERE
    return np.exp(-0.5 * a ** 2) * np.cosh(a)

```

```python
def phi_2(a):
    return -a ** 3

def p_2(a):
    # YOUR CODE HERE
    return np.exp(-0.25 * a ** 4)

```

```python
import math

def phi_3(a):
    return -6 * a / (a ** 2 + 5)

def p_3(a):
    # YOUR CODE HERE
    # = exp( - integral(-6*a / (a**2 + 5)) ~> exp(-3 ln (a^2+5)) ~> (a^2+5)^(-3)
    return np.power((np.power(a, 2) + 5), -3)
    # return (a ** 2 + 5) ** (-3)

```

```python
activation_functions = [phi_0, phi_1, phi_2, phi_3]
priors = [p_0, p_1, p_2, p_3]

a = np.linspace(-5, 5, 1000)
for prior in priors:
    assert prior(a).shape == (1000,), "Wrong output shape"

```

Plot the activation functions and the corresponding prior distributions, from $a = -5$ to $5$ (hint: use the lists defined in the cell above). Compare the shape of the priors to the histogram you plotted in the last question.

```python
### 1.3 Implicit priors (continued)
a = np.linspace(-5, 5, 1000)

for index, funct in enumerate(activation_functions):
    plt.subplot(1, len(activation_functions), index + 1)
    plt.plot(a, funct(a))
    plt.suptitle('Activation functions')
```

```python
for index, funct in enumerate(priors):
    plt.subplot(1, len(priors), index + 1)
    plt.plot(a, funct(a))
    plt.suptitle('Priors')
```

### Answer to: Compare the shape of the priors to the histogram you plotted in the last question.
We can see that all priors somewhat resemble the shape of a Gaussian distributions. Therefore, they are all similar to the shape of the histograms of the measurements, as we have identified them to be close to a Gaussian distribution as well. However, we notice that priors 1 and 2 show a broader peak, which is also visible in the aforementioned histograms. Therefore, we conclude, that they resemble the given measurements the best.

### 1.4 Whitening (15 points)
Some ICA algorithms can only learn from whitened data. Write a method `whiten(X)` that takes a $M \times T$ data matrix $\mathbf{X}$ (where $M$ is the dimensionality and $T$ the number of examples) and returns a whitened matrix. If you forgot what whitening is or how to compute it, various good sources are available online, such as http://courses.media.mit.edu/2010fall/mas622j/whiten.pdf. Your function should also center the data before whitening.

```python
### 1.4 Whitening
def whiten(X):
    # normalize the mean to zero
    X -= np.mean(X, axis=1, keepdims=True)
    # calculate the covariance matrix
    Phi = np.cov(X)
    # calculate eigenvalues and eigenvalues
    Lambda, Phi = np.linalg.eigh(Phi)
    # set Lambda to have the eigenvalues on the diagonal
    Lambda = np.diag(Lambda)
    # calculate A
    A = Phi.dot(np.sqrt(np.linalg.inv(Lambda)))
    # get whitend data
    w = A.T.dot(X)
    return w

```

```python
### Test your function
Xw = whiten(X)
assert Xw.shape == (num_sources, signal_length), "The shape of your mixed signals is incorrect"

```

### 1.5 Interpret results of whitening (10 points)
Make 3 figures, one for the sources, one for measurements and one for the whitened measurements. In each figure, make $5 \times 5$ subplots with scatter plots for each pair of signals. Each axis represents a signal and each time-instance is plotted as a dot in this space. You can use the `plt.scatter()` function. Describe what you see.

Now compute and visualize the covariance matrix of the sources, the measurements and the whitened measurements. You can visualize each covariance matrix using this code:
```python
# Dummy covariance matrix C;
C = np.eye(5)
ax = imshow(C, cmap='gray', interpolation='nearest')
```

```python
### 1.5 Interpret results of whitening
def scatter_plot(signal):
    fig, axes = plt.subplots(len(signal), len(signal))

    for i, signal_1 in enumerate(signal):
        for j, signal_2 in enumerate(signal):

            axes[i][j].scatter(signal_1, signal_2, s = 0.2)

```

```python
# plot signals

scatter_plot(S)
```

```python
# plot measure
scatter_plot(X)
```

```python
# plot whitend measure
scatter_plot(Xw)
```

```python
# plot the covariance matices

for index, conv in enumerate([np.cov(S), np.cov(X), np.cov(Xw)]):
    plt.subplot(1, 3, index + 1)
    imshow(conv, cmap='gray', interpolation='nearest')
    #print(conv)
```

Are the signals independent after whitening?

YOUR ANSWER HERE

No, they are decorrelated but not necessarily independent.

### 1.6 Covariance (5 points)
Explain what a covariant algorithm is.

YOUR ANSWER HERE

A covariant algorithm is a consistent algorithm that should give the same results independent of the units in which quantities are measured (Knuth, 1968).

If we look at the steepest descent rule, for example, we see that it is a non-covariant algorithm. Its update rule consists of a term of dimensions $[w_i]$ that is set equal to a term of dimensions $1/[w_i]$. This leads to the fact that this non-covariant algorithm needs a very large number of iterations to achieve convergence.
In order to avoid this, we can adjust the update rule by multiplying the right-hand side with a matrix M of dimensions $[w_i, 1/w_i]$ in order to get the same dimensionality on both sides of the equation. This gives a covariant algorithm. M might be found by looking at the curvature of the objective function.

### 1.7 Independent Component Analysis (25 points)
Implement the covariant ICA algorithm as described in MacKay. Write a function `ICA(X, activation_function, learning_rate)`, that returns the demixing matrix $\mathbf{W}$. The input `activation_function` should accept a function such as `lambda a: -tanh(a)`. Update the gradient in batch mode, averaging the gradients over the whole dataset for each update. Make it efficient, so use matrix operations instead of loops where possible (loops are slow in interpreted languages such as python and matlab, whereas matrix operations are internally computed using fast C code). Experiment with the learning rate and the initialization of $\mathbf{W}$. Your algorithm should be able to converge (i.e. `np.linalg.norm(grad) < 1e-5`) in within 10000 steps.

```python
def ICA(X, activation_function, learning_rate=0.01, num_iterations = 10000):
    # YOUR CODE HERE
    # initialize W randomly
    W = random_nonsingular_matrix(d=len(X))

    # only loop over iterations
    for _ in range(num_iterations):
        # 1) Put x through a linear mapping:
        a = W @ X
        # 2) Put a through a nonlinear map:
        z = activation_function(a)
        # 3) Put a back through W:
        x_ = W.T @ a
        # calculate gradient
        grad = W + z @ x_.T / X.shape[1]
        # check for convergence
        if np.linalg.norm(grad) < 1e-5:
            #print(np.linalg.norm(grad))
            break
        # 4) Adjust the weights in accordance with
        W += learning_rate * grad
    return W

W_est = ICA(Xw, phi_3)  # Compare with ICA(X)

## Testing different learning rates, we found that the value 0.01 works the best.
## With higher values we experience numerical overflow in the square of phi_3.

## We found that initializing W with a uniform distribution (using numpy.random.rand and learning rate 0.01)
## works better than when using a Gaussian distribution (using numpy.random.randn and learning rate 0.001).
```

```python
# We will test your function so make sure it runs with only X and phi as input, and returns only W
# Also it should converge for all activation functions

W_estimates = [ICA(Xw, activation_function=phi) for phi in activation_functions]
assert all([W_est.shape == (num_sources, num_sources) for W_est in W_estimates])

```

### 1.8 Experiments  (5 points)
Run ICA on the provided signals using each activation function $\phi_0, \ldots, \phi_3$ (or reuse `W_estimates`). Use the found demixing matrix $\mathbf{W}$ to reconstruct the signals and plot the retreived signals for each choice of activation function.

```python
# 1.8 Experiments
def plot_results(W_estimates, X):
  for index, W in enumerate(W_estimates):
      plot_signals(W @ X, 'activation: Phi_{}'.format(index))

plot_results(W_estimates, Xw)
```

```python
import time

def does_whitening_make_a_difference():
    # Does it make a difference (in terms of speed of convergence)
    # if you whiten your data before running ICA?

    # Return True or False
    # YOUR CODE HERE
    t0_w = time.time()
    ## increasing the number of possible iterations to enable convergence and get meaningful results
    W_estimates_w = [ICA(Xw, activation_function=phi,num_iterations=100000) for phi in activation_functions]
    t1_w = time.time()

    t0 = time.time()
    W_estimates = [ICA(X, activation_function=phi, num_iterations=100000) for phi in activation_functions]
    t1 = time.time()

    time_whitened = t1_w - t0_w
    time_non_whitened = t1 - t0
    print('time (sec) using whitened data: ', time_whitened)
    print('time (sec) using the original data: ', time_non_whitened)
    print('')

    if time_non_whitened / time_whitened > 1.2:
      print('Whitening does make a difference.')
      return True
    else:
      print('Whitening does not make a difference.')
      return False

```

```python
assert type(does_whitening_make_a_difference()) == bool
```

### 1.9 Audio demixing (10 points)
The 'cocktail party effect' refers to the ability humans have to attend to one speaker in a noisy room. We will now use ICA to solve a similar but somewhat idealized version of this problem. The code below loads 5 sound files and plots them.

Use a random non-singular mixing matrix to mix the 5 sound files. You can listen to the results in your browser using `play_signals`, or save them to disk if this does not work for you. Plot histograms of the mixed audio and use your ICA implementation to de-mix these and reproduce the original source signals. As in the previous exercise, try each of the activation functions.

Keep in mind that this problem is easier than the real cocktail party problem, because in real life there are often more sources than measurements (we have only two ears!), and the number of sources is unknown and variable. Also, mixing is not instantaneous in real life, because the sound from one source arrives at each ear at a different point in time. If you have time left, you can think of ways to deal with these issues.

```python
import scipy.io.wavfile
from IPython.display import Audio, display, Markdown

# Save mixtures to disk, so you can listen to them in your audio player
def save_wav(data, out_file, rate):
    scaled = np.int16(data / np.max(np.abs(data)) * 32767)
    scipy.io.wavfile.write(out_file, rate, scaled)

# Or play them in your browser
def play_signals(S, sample_rate, title="Signals"):
    display(Markdown(title))
    for signal in S:
        display(Audio(signal, rate=sample_rate))
```

```python
# Load audio sources
source_files = ['beet.wav', 'beet9.wav', 'beet92.wav', 'mike.wav', 'street.wav']

wav_data = []
sample_rate = None
for f in source_files:
    sr, data = scipy.io.wavfile.read(f, mmap=False)
    if sample_rate is None:
        sample_rate = sr
    else:
        assert(sample_rate == sr)
    wav_data.append(data[:190000])  # cut off the last part so that all signals have same length

# Create source and measurement data
S_audio = np.c_[wav_data]
plot_signals(S_audio)
# play_signals(S_audio, sample_rate)
```

```python
### 1.9 Audio demixing
X_mix = make_mixtures(S_audio, random_nonsingular_matrix(d=len(S)))

# whiten data
X_mix_w = whiten(X_mix)

# find W estimates with different phi and plot results
W_estimates_audio = [ICA(X_mix_w, activation_function=phi) for phi in activation_functions]
```

```python
plot_signals(X_mix, 'Generated mixture')
```

```python
plot_histograms(X_mix, 'Histogram of mixed audio signals')

plot_results(W_estimates_audio, X_mix_w)

reconstructed_S = W_estimates_audio @ X_mix_w
play_signals(reconstructed_S[3], sample_rate)
```

Report your results. Using which activation functions ICA recovers the sources?

YOUR ANSWER HERE
Again, we can see that the mixed signals take a form that is close to that of a Gaussian distribution.
Both by examining the plots above and listening to the results, we can conclude that the activation functions 0 and 3 manage to separate the mixed signals into the original signals the best.

### 1.10 Excess Kurtosis (15 points)
The (excess) kurtosis is a measure of 'peakedness' of a distribution. It is defined as
$$
\verb+Kurt+[X] = \frac{\mu_4}{\sigma^4} - 3 = \frac{\operatorname{E}[(X-{\mu})^4]}{(\operatorname{E}[(X-{\mu})^2])^2} - 3
$$
Here, $\mu_4$ is known as the fourth moment about the mean, and $\sigma$ is the standard deviation.
The '-3' term is introduced so that a Gaussian random variable has 0 excess kurtosis.
We will now try to understand the performance of the various activation functions by considering the kurtosis of the corresponding priors, and comparing those to the empirical kurtosis of our data.

#### 1.10.1 (10 points)
First, compute analytically the kurtosis of the four priors that you derived from the activation functions before. You may find it helpful to use an online service such as [Wolfram Alpha](https://www.wolframalpha.com/) or [Integral Calculator](https://www.integral-calculator.com/) to (help you) evaluate the required integrals. Give your answer as both an exact expression as well as a numerical approximation (for example $\frac{\pi}{2} \approx 1.571$).

To calculate the moments we first need the mean and the prior. For the mean we see that all priors are zero-centered and therefore we have a mean of 0. Using the results obtained in the exercises above we can get the kurtosis as follows:

Case 1:
\begin{align}
p(a) & = \frac{1}{\pi cosh(a)}\\
\mathbb{E}[X^4] &= \int \frac{1}{\pi cosh(a)}x^4 dx\\
& = \frac{5\pi^4}{16}\\
\mathbb{E}[X^2]^2 &= \left[\int \frac{1}{\pi cosh(a)}x^2  dx\right]^2 \\
& = \frac{\pi^4}{16}\\
Kurt[X] &= \frac{\frac{5\pi^4}{16}}{\frac{\pi^4}{16}}-3\\
& = 2
\end{align}

Case 2:
\begin{align}
p(a) & = \frac{1}{\sqrt{2\pi e}}exp\left(\frac{-a^2}{2}\right)cosh(a)\\
\mathbb{E}[X^4] &= \int \frac{1}{\sqrt{2\pi e}}exp\left(\frac{-a^2}{2}\right)cosh(a)  x^4 dx\\
& = 10 \frac{\sqrt{2\pi e}}{\sqrt{2\pi e}}\\
\mathbb{E}[X^2]^2 &= \left[\int \frac{1}{\sqrt{2\pi e}}exp\left(\frac{-a^2}{2}\right)cosh(a)  x^2 dx\right]^2 \\
& = \left(2\frac{\sqrt{2\pi e}}{\sqrt{2\pi e}}\right)^2\\
Kurt[X] &= \frac{ 10 \frac{\sqrt{2\pi e}}{\sqrt{2\pi e}}}{\left(2\frac{\sqrt{2\pi e}}{\sqrt{2\pi e}}\right)^2}-3\\
&=-0.5
\end{align}

Case 3:
\begin{align}
p(a) & = \frac{\sqrt{2}}{\Gamma(0.25)}exp\left(-\frac{x^4}{4}\right)\\
\mathbb{E}[X^4] &= \int \frac{\sqrt{2}}{\Gamma(0.25)}exp\left(-\frac{x^4}{4}\right) x^4 dx\\
& =  \frac{\sqrt{2}}{\Gamma(0.25)}\frac{\Gamma(0.25)}{\sqrt{2}}\\
\mathbb{E}[X^2]^2 &= \left[\int \frac{\sqrt{2}}{\Gamma(0.25)}exp\left(-\frac{x^4}{4}\right)  x^2 dx\right]^2 \\
& =  \frac{2}{\Gamma^2(0.25)}\left( \sqrt{2} \Gamma(3/4)\right)^2\\
Kurt[X] &= \frac{\frac{\sqrt{2}}{\Gamma(0.25)}\frac{\Gamma(0.25)}{\sqrt{2}}}{\frac{2}{\Gamma^2(0.25)}\left( \sqrt{2} \Gamma(3/4)\right)^2} -3\\
&=\frac{1}{\frac{4\Gamma^2(3/4)}{\Gamma^2(1/4)}}-3\\
&=\frac{\Gamma^2(1/4)}{4\Gamma^2(3/4)}-3\\
&\approx -0.811
\end{align}

Case 4:
\begin{align}
p(a) & = \frac{3\pi}{200\sqrt{5}}(a^2+5)^{-3}\\
\mathbb{E}[X^4] &= \int \frac{3\pi}{200\sqrt{5}}(x^2+5)^{-3} x^4 dx\\
\mathbb{E}[X^2]^2 &= \left[\int \frac{3\pi}{200\sqrt{5}}(a^2+5)^{-3}  x^2 dx\right]^2 \\
Kurt[X] &= \frac{\int \frac{3\pi}{200\sqrt{5}}(x^2+5)^{-3} x^4 dx}{\left[\int \frac{3\pi}{200\sqrt{5}}(a^2+5)^{-3}  x^2 dx\right]^2}-3\\
&= 6
\end{align}

```python
### Include your answer here (you can use math.gamma if needed)
def get_kurtosis():
    # Return a list with 4 numbers / expressions
    # return [0, 0, 0, 0]

    k1 = [
        2,
        -0.5,
        np.power(math.gamma(1/4), 2)/(np.power(math.gamma(3/4), 2) * 4) - 3,
        6
    ]

    return k1
```

```python
# Let's check
kurtosis = get_kurtosis()
print (kurtosis)
assert len(kurtosis) == 4
```

#### 1.10.2 (5 points)
Now use the `scipy.stats.kurtosis` function, with the `fisher` option set to `True`, to compute the empirical kurtosis of the dummy signals and the real audio signals.

Can you use this data to explain the performance of the various activation functions on the synthetic and real data?

```python
### 1.10.2 Excess Kurtosis
# YOUR CODE HERE

import scipy.stats

def get_kurtosis_empirical():
    # theoretical benchmark

    #k1 = get_kurtosis()
    #print('Analytical Kurtosis: {}\n\n'.format(k1))

    k2 = scipy.stats.kurtosis(S_audio, axis=1, fisher=True)
    k3 = scipy.stats.kurtosis(reconstructed_S, axis=2, fisher=True)
    k4 = scipy.stats.kurtosis(X, axis=1, fisher=True)
    reconstructed_X = W_estimates @ Xw
    k5 = scipy.stats.kurtosis(reconstructed_X, axis=2, fisher=True)

    l = ['original audio', 'reconstructed audio (activation functions 0 to 3)', 'original synthetic data', 'reconstructed synthetic data (activation functions 0 to 3)']

    for index, k in enumerate([k2, k3, k4, k5]):
        print('Kurtosis of {}:\n {}\n'.format(l[index], k))

get_kurtosis_empirical()

```

YOUR ANSWER HERE

Yes, we can.

In general, the kurtosis is a measure for outlayers: values around the mean are not relevant and do not contribute to the overall value of the kurtosis. Therefore, it is a measure of outlayers / the tails.

From the empirical kurtosis values, we can draw the same conclusions about the performances of the different activation functions as seen above. For the synthetic data, we see that the activation functions 0 and 3 both return kurtosis values that are all close to zero. The activation functions 1 and 2, that also showed the best separation of the mixed signals as shown in the plots above, manage to get kurtosis value similar to the original data. Note, that the values are scrambled, but it is possible to find all the corresponding kurtosis values in the reconstructed data.

Comparing the kurtosis values for the different activation functions when applying them to the actual audio files, we find that again they match our previous observation that the function 0 and 3 manage to reconstruct the original signal the best. Even though they are scrambled, activation function 0 and 3 show very similar kurtosis values as the original audio, while functions 1 and 2 show numbers close to zero.


text:
# Tutorial 2: Differential Equations
**Week 0, Day 4: Calculus**

**By Neuromatch Academy**

__Content creators:__ John S Butler, Arvind Kumar with help from Rebecca Brady

__Content reviewers:__  Swapnil Kumar, Sirisha Sripada, Matthew McCann, Tessy Tom

__Production editors:__ Matthew McCann, Ella Batty

**Our 2021 Sponsors, including Presenting Sponsor Facebook Reality Labs**

---
# Tutorial Objectives

A great deal of neuroscience can be modelled using differential equations, from gating channels to single neurons to a network of neurons to blood flow, to behaviour. A simple way to think about differential equations is they are equations that describe how something changes.

The most famous of these in neuroscience is the Nobel Prize winning Hodgkin Huxley equation, which describes a neuron by modelling the gating of each axon. But we will not start there; we will start a few steps back.

Differential Equations are mathematical equations that describe how something like population or a neuron changes over time. The reason why differential equations are so useful is they can generalise a process such that one equation can be used to describe many different outcomes.
The general form of a first order differential equation is:

\begin{align*}
\frac{d}{dt}y(t)&=f(t,y(t))\\
\end{align*}

which can be read as "the change in a process $y$ over time $t$ is a function $f$ of time $t$ and itself $y$". This might initially seem like a paradox as you are using a process $y$ you want to know about to describe itself, a bit like the MC Escher drawing of two hands painting [each other](https://en.wikipedia.org/wiki/Drawing_Hands). But that is the beauty of mathematics - this can be solved some of time, and when it cannot be solved exactly we can use numerical methods to estimate the answer (as we will see in the next tutorial).

In this tutorial, we will see how __differential equations are motivated by observations of physical responses.__ We will break down the population differential equation, then the integrate and fire model, which leads nicely into raster plots and frequency-current curves to rate models.

**Steps:**
- Get an intuitive understanding of a linear population differential equation (humans, not neurons)
- Visualize the relationship between the change in population and the population
- Breakdown the Leaky Integrate and Fire (LIF) differential equation
- Code the exact solution of an LIF for a constant input
- Visualize and listen to the response of the LIF for different inputs

##  Video 1: Why do we care about differential equations?

```python
# @title Video 1: Why do we care about differential equations?
from ipywidgets import widgets

out2 = widgets.Output()
with out2:
  from IPython.display import IFrame
  class BiliVideo(IFrame):
      def __init__(self, id, page=1, width=400, height=300, **kwargs):
          self.id=id
          src = 'https://player.bilibili.com/player.html?bvid={0}&page={1}'.format(id, page)
          super(BiliVideo, self).__init__(src, width, height, **kwargs)

  video = BiliVideo(id="BV1v64y197bW", width=730, height=410, fs=1)
  print('Video available at https://www.bilibili.com/video/{0}'.format(video.id))
  display(video)

out1 = widgets.Output()
with out1:
  from IPython.display import YouTubeVideo
  video = YouTubeVideo(id="LhX-mUd8lPo", width=730, height=410, fs=1, rel=0)
  print('Video available at https://youtube.com/watch?v=' + video.id)
  display(video)

out = widgets.Tab([out1, out2])
out.set_title(0, 'Youtube')
out.set_title(1, 'Bilibili')

display(out)
```

    Tab(children=(Output(), Output()), _titles={'0': 'Youtube', '1': 'Bilibili'})

---
# Setup

```python
# Imports

import numpy as np
import matplotlib.pyplot as plt
```

##  Figure Settings

```python
# @title Figure Settings
import IPython.display as ipd
from matplotlib import gridspec

import ipywidgets as widgets  # interactive display
%config InlineBackend.figure_format = 'retina'

# use NMA plot style
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")
my_layout = widgets.Layout()
```

##  Plotting Functions

```python
# @title Plotting Functions

def plot_dPdt(alpha=.3):
  """ Plots change in population over time
    Args:
      alpha: Birth Rate
    Returns:
      A figure two panel figure
      left panel: change in population as a function of population
      right panel: membrane potential as a function of time
  """

  with plt.xkcd():
    time=np.arange(0, 10 ,0.01)
    fig = plt.figure(figsize=(12,4))
    gs = gridspec.GridSpec(1, 2)

    ## dpdt as a fucntion of p
    plt.subplot(gs[0])
    plt.plot(np.exp(alpha*time), alpha*np.exp(alpha*time))
    plt.xlabel(r'Population $p(t)$ (millions)')
    plt.ylabel(r'$\frac{d}{dt}p(t)=\alpha p(t)$')

    ## p exact solution
    plt.subplot(gs[1])
    plt.plot(time, np.exp(alpha*time))
    plt.ylabel(r'Population $p(t)$ (millions)')
    plt.xlabel('time (years)')
    plt.show()

def plot_V_no_input(V_reset=-75):
  """
    Args:
      V_reset: Reset Potential
    Returns:
      A figure two panel figure
      left panel: change in membrane potential as a function of membrane potential
      right panel: membrane potential as a function of time
  """
  E_L=-75
  tau_m=10
  t=np.arange(0,100,0.01)
  V= E_L+(V_reset-E_L)*np.exp(-(t)/tau_m)
  V_range=np.arange(-90,0,1)
  dVdt=-(V_range-E_L)/tau_m

  with plt.xkcd():
    time=np.arange(0, 10, 0.01)
    fig = plt.figure(figsize=(12, 4))
    gs = gridspec.GridSpec(1, 2)

    plt.subplot(gs[0])
    plt.plot(V_range,dVdt)
    plt.hlines(0,min(V_range),max(V_range), colors='black', linestyles='dashed')
    plt.vlines(-75, min(dVdt), max(dVdt), colors='black', linestyles='dashed')
    plt.plot(V_reset,-(V_reset - E_L)/tau_m, 'o', label=r'$V_{reset}$')
    plt.text(-50, 1, 'Positive')
    plt.text(-50, -2, 'Negative')
    plt.text(E_L - 1, max(dVdt), r'$E_L$')
    plt.legend()
    plt.xlabel('Membrane Potential V (mV)')
    plt.ylabel(r'$\frac{dV}{dt}=\frac{-(V(t)-E_L)}{\tau_m}$')

    plt.subplot(gs[1])
    plt.plot(t,V)
    plt.plot(t[0],V_reset,'o')
    plt.ylabel(r'Membrane Potential $V(t)$ (mV)')
    plt.xlabel('time (ms)')
    plt.ylim([-95, -60])

    plt.show()

## LIF PLOT
def plot_IF(t, V,I,Spike_time):
  """
    Args:
      t  : time
      V  : membrane Voltage
      I  : Input
      Spike_time : Spike_times
    Returns:
      figure with three panels
      top panel: Input as a function of time
      middle panel: membrane potential as a function of time
      bottom panel: Raster plot
  """

  with plt.xkcd():
    fig = plt.figure(figsize=(12, 4))
    gs = gridspec.GridSpec(3, 1,  height_ratios=[1, 4, 1])

    # PLOT OF INPUT
    plt.subplot(gs[0])
    plt.ylabel(r'$I_e(nA)$')
    plt.yticks(rotation=45)
    plt.hlines(I,min(t),max(t),'g')
    plt.ylim((2, 4))
    plt.xlim((-50, 1000))

    # PLOT OF ACTIVITY
    plt.subplot(gs[1])
    plt.plot(t,V)
    plt.xlim((-50, 1000))
    plt.ylabel(r'$V(t)$(mV)')

    # PLOT OF SPIKES
    plt.subplot(gs[2])
    plt.ylabel(r'Spike')
    plt.yticks([])
    plt.scatter(Spike_time, 1 * np.ones(len(Spike_time)), color="grey", marker=".")
    plt.xlim((-50, 1000))
    plt.xlabel('time(ms)')
    plt.show()

##    Plotting the differential Equation
def plot_dVdt(I=0):
  """
    Args:
      I  : Input Current
    Returns:
      figure of change in membrane potential as a function of membrane potential
  """

  with plt.xkcd():
    E_L = -75
    tau_m = 10
    V = np.arange(-85, 0, 1)
    g_L = 10.
    fig = plt.figure(figsize=(6, 4))

    plt.plot(V,(-(V-E_L) + I*10) / tau_m)
    plt.hlines(0, min(V), max(V), colors='black', linestyles='dashed')
    plt.xlabel('V (mV)')
    plt.ylabel(r'$\frac{dV}{dt}$')
    plt.show()
```

##  Helper Functions

```python
# @title Helper Functions

## EXACT SOLUTION OF LIF
def Exact_Integrate_and_Fire(I,t):
  """
    Args:
      I  : Input Current
      t : time
    Returns:
      Spike : Spike Count
      Spike_time : Spike time
      V_exact : Exact membrane potential
  """

  Spike = 0
  tau_m = 10
  R = 10
  t_isi = 0
  V_reset = E_L = -75
  V_exact = V_reset * np.ones(len(t))
  V_th = -50
  Spike_time = []

  for i in range(0, len(t)):

    V_exact[i] = E_L + R*I + (V_reset - E_L - R*I) * np.exp(-(t[i]-t_isi)/tau_m)

    # Threshold Reset
    if V_exact[i] > V_th:
        V_exact[i-1] = 0
        V_exact[i] = V_reset
        t_isi = t[i]
        Spike = Spike+1
        Spike_time = np.append(Spike_time, t[i])

  return Spike, Spike_time, V_exact
```

---
# Section 1: Population differential equation

##  Video 2: Population differential equation

```python
# @title Video 2: Population differential equation
from ipywidgets import widgets

out2 = widgets.Output()
with out2:
  from IPython.display import IFrame
  class BiliVideo(IFrame):
      def __init__(self, id, page=1, width=400, height=300, **kwargs):
          self.id=id
          src = 'https://player.bilibili.com/player.html?bvid={0}&page={1}'.format(id, page)
          super(BiliVideo, self).__init__(src, width, height, **kwargs)

  video = BiliVideo(id="BV1pg41137CU", width=730, height=410, fs=1)
  print('Video available at https://www.bilibili.com/video/{0}'.format(video.id))
  display(video)

out1 = widgets.Output()
with out1:
  from IPython.display import YouTubeVideo
  video = YouTubeVideo(id="czgGyoUsRoQ", width=730, height=410, fs=1, rel=0)
  print('Video available at https://youtube.com/watch?v=' + video.id)
  display(video)

out = widgets.Tab([out1, out2])
out.set_title(0, 'Youtube')
out.set_title(1, 'Bilibili')

display(out)
```

    Tab(children=(Output(), Output()), _titles={'0': 'Youtube', '1': 'Bilibili'})

This video covers our first example of a differential equation: a differential equation which models the change in population.

 Click here for text recap of video

To get an intuitive feel of a differential equations, we will start with a population differential equation, which models the change in population [1], that is human population not neurons, we will get to neurons later. Mathematically it is written like:
\begin{align*}
\\
\frac{d}{dt}\,p(t) &= \alpha p(t),\\
\end{align*}

where $p(t)$ is the population of the world and $\alpha$ is a parameter representing birth rate.

Another way of thinking about the models is that the equation
\begin{align*}
\\
\frac{d}{dt}\,p(t) &= \alpha p(t),\\
\text{can be written as:}\\
\text{"Change in Population"} &= \text{ "Birth rate times Current population."}
\end{align*}

The equation is saying something reasonable maybe not the perfect model but a good start.

### Think! 1.1: Interpretating the behavior of a linear population equation
Using the plot below of change of population $\frac{d}{dt} p(t) $ as a function of population $p(t)$ with birth-rate $\alpha=0.3$, discuss the following questions:
1. Why is the population differential equation known as a linear differential equation?

2. How does population size affect the rate of change of the population?

```python
# @markdown Execute the code to plot the rate of change of population as a function of population
p = np.arange(0, 100, 0.1)

with plt.xkcd():

  dpdt = 0.3*p
  fig = plt.figure(figsize=(6, 4))
  plt.plot(p, dpdt)
  plt.xlabel(r'Population $p(t)$ (millions)')
  plt.ylabel(r'$\frac{d}{dt}p(t)=\alpha p(t)$')
  plt.show()
```

[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D4_Calculus/solutions/W0D4_Tutorial2_Solution_d65c50d6.py)

1. The plot of $\frac{dp}{dt}$ is a line, which is why the differential
  equation is known as a linear differential equation.

2. As the population increases, the change of population increases. A
  population of 20 has a change of 6 while a population of 100 has a change of
  30. This makes sense - the larger the population the larger the change.

## Section 1.1: Exact solution of the population equation

### Section 1.1.1: Initial condition
The linear population differential equation is known as an initial value differential equation because we need an initial population value to solve it. Here we will set our initial population to:

\begin{align*}
&p(0)=1.\\
\end{align*}

Different initial conditions will lead to different answers, but they will not change the differential equation. This is one of the strengths of a differential equation.

###  Section 1.1.2: Exact Solution
To calculate the exact solution of a differential equation, we must integrate both sides. Instead of numerical integration (as you delved into in the last tutorial), we will first try to solve the differential equations using analytical integration. As with derivatives, we can find analytical integrals of simple equations by consulting [a list](https://en.wikipedia.org/wiki/Lists_of_integrals). We can then get intregrals for more complex equations using some mathematical tricks - the harder the equation the more obscure the trick.

Given the linear population equation
\begin{align*}
\frac{d}{dt}\,p(t) &= \alpha p(t),\\\\
p(0)=P_0,\\
\end{align*}
has the exact solution:
\begin{align*}
p(t)&=P_0e^{\alpha t}.\\
\end{align*}

The exact solution written in words is:

\begin{align*}
\text{"Population"}&=\text{"grows/declines exponentially as a function time and birth rate"}.\\
\end{align*}

Most differential equations do not have a known exact solution, so in the next tutorial on numerical methods we will show how the solution can be estimated.

A small aside: a good deal of progress in mathematics was due to mathematicians writing taunting letters to each other saying they had a trick that could solve something better than everyone else. So do not worry too much about the tricks.

####  Example Exact Solution of the Population Equation
Let's consider the population differential equation with a birth rate $\alpha=0.3$:

\begin{align*}
\frac{d}{dt}\,p(t) = 0.3 p(t),\\
\text{with the initial condition}\\
p(0)=1.\\
\end{align*}

It has an exact solution
\begin{align*}
\\
p(t)=e^{0.3 t}.
\end{align*}

```python
# @markdown Execute code to plot the exact solution
t = np.arange(0, 10, 0.1) # Time from 0 to 10 years in 0.1 steps

with plt.xkcd():

  p = np.exp(0.3 * t)

  fig = plt.figure(figsize=(6, 4))
  plt.plot(t, p)
  plt.ylabel('Population (millions)')
  plt.xlabel('time (years)')
  plt.show()
```

## Section 1.2: Parameters of the differential equation
One of the goals when designing a differential equation is to make it generalisable. Which means that the differential equation will give reasonable solutions for different countries with different birth rates $\alpha$.

###  Interactive Demo 1.2:  Interactive Parameter Change
Play with the widget to see the relationship between $\alpha$ and the population differential equation as a function of population (left-hand side), and the population solution as a function of time (right-hand side). Pay close attention to the transition point from positive to negative.

How do changing parameters of the population equation affect the outcome?

1. What happens when $\alpha < 0$?
2. What happens when $\alpha > 0$?
3. What happens when $\alpha = 0$?

```python
# @markdown Make sure you execute this cell to enable the widget!
my_layout.width = '450px'
@widgets.interact(
    alpha=widgets.FloatSlider(.3, min=-1., max=1., step=.1, layout=my_layout)
)
def Pop_widget(alpha):
  plot_dPdt(alpha=alpha)
  plt.show()
```

    interactive(children=(FloatSlider(value=0.3, description='alpha', layout=Layout(width='450px'), max=1.0, min=-…

[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D4_Calculus/solutions/W0D4_Tutorial2_Solution_7cf110c5.py)

1. Negative values of alpha result in an exponential decrease to 0 a stable solution.

2. Positive Values of alpha in an exponential increases to infinity.

3. Alpha equal to 0 is a unique point known as an equilibrium point when the
  dp/dt=0 and there is no change in population. This is known as a stable point.

The population differential equation is an over-simplification and has some very obvious limitations:
1. Population growth is not exponential as there are limited number of resources so the population will level out at some point.
2. It does not include any external factors on the populations like weather, predators and preys.

These kind of limitations can be addressed by extending the model.

While it might not seem that the population equation has direct relevance to neuroscience a similar equation is used to describe the accumulation of evidence for decision making known as the Drift Diffusion Model which you will see in more detail in the Linear System day in Neuromatch (W2D2).

Another differential equation that is similar to the population equation is the Leaky Integrate and Fire model which you have seen in the python pre-course materials on W0D1 and W0D2, it will turn up later in Neuromatch as well. Below we will delve in the motivation of the differential equation.

---
# Section 2: The leaky integrate and fire model

##  Video 3: The leaky integrate and fire model

```python
# @title Video 3: The leaky integrate and fire model
from ipywidgets import widgets

out2 = widgets.Output()
with out2:
  from IPython.display import IFrame
  class BiliVideo(IFrame):
      def __init__(self, id, page=1, width=400, height=300, **kwargs):
          self.id=id
          src = 'https://player.bilibili.com/player.html?bvid={0}&page={1}'.format(id, page)
          super(BiliVideo, self).__init__(src, width, height, **kwargs)

  video = BiliVideo(id="BV1rb4y1C79n", width=730, height=410, fs=1)
  print('Video available at https://www.bilibili.com/video/{0}'.format(video.id))
  display(video)

out1 = widgets.Output()
with out1:
  from IPython.display import YouTubeVideo
  video = YouTubeVideo(id="ZfWO6MLCa1s", width=730, height=410, fs=1, rel=0)
  print('Video available at https://youtube.com/watch?v=' + video.id)
  display(video)

out = widgets.Tab([out1, out2])
out.set_title(0, 'Youtube')
out.set_title(1, 'Bilibili')

display(out)
```

    Tab(children=(Output(), Output()), _titles={'0': 'Youtube', '1': 'Bilibili'})

This video covers the Leaky Integrate and Fire model (a linear differential equation which describes the membrane potential of a single neuron).

 Click here for text recap of full LIF equation from video

The Leaky Integrate and Fire Model is a linear differential equation that describes the membrane potential ($V$) of a single neuron which was proposed by Louis Édouard Lapicque in 1907 [2].

The subthreshold membrane potential dynamics of a LIF neuron is described by
\begin{align}
\tau_m\frac{dV}{dt} = -(V-E_L) + R_mI\,
\end{align}

where $\tau_m$ is the time constant, $V$ is the membrane potential,  $E_L$ is the resting potential, $R_m$ is membrane resistance, and $I$ is the external input current.

In the next few sections, we will break down the full LIF equation and then build it back up to get an intuitive feel of the different facets of the differential equation.

## Section 2.1: LIF without input

As seen in the video, we will first model an LIF neuron without input, which results in the equation:

\begin{align}
\frac{dV}{dt} &= \frac{-(V-E_L)}{\tau_m}.\\
\end{align}

where $\tau_m$ is the time constant, $V$ is the membrane potential, and $E_L$ is the resting potential.

 Click here for further details (from video)

First, we will simplify the equation by removing the input, which gives the equation

\begin{align}
\tau_m\frac{dV}{dt} &= -V+E_L,\\
\end{align}

which can be written in words as:

\begin{align}
\begin{matrix}\text{"Time constant multiplied by the} \\ \text{change in membrane potential"}\end{matrix}&=\begin{matrix}\text{"Minus Current} \\ \text{membrane potential"} \end{matrix}+
\begin{matrix}\text{"resting potential"}\end{matrix}.\\
\end{align}

The equation can be re-arranged to look even more like the population equation:

\begin{align}
\frac{dV}{dt} &= \frac{-(V-E_L)}{\tau_m}.\\
\end{align}

### Think! 2.1: Effect on membrane potential $V$ on the LIF model

The plot the below shows the change in membrane potential $\frac{dV}{dt}$ as a function of membrane potential $V$ with the parameters set as:
* `E_L = -75`
* `V_reset = -50`
* `tau_m = 10.`

1. What is the effect on $\frac{dV}{dt}$ when $V>-75$ mV?
2. What is the effect on $\frac{dV}{dt}$ when $V<-75$ mV
3. What is the effect on $\frac{dV}{dt}$ when $V=-75$ mV?

```python
# @markdown Make sure you execute this cell to plot the relationship between dV/dt and V
# Parameter definition
E_L = -75
tau_m = 10

# Range of Values of V
V = np.arange(-90, 0, 1)
dV = -(V - E_L) / tau_m

with plt.xkcd():

  fig = plt.figure(figsize=(6, 4))
  plt.plot(V, dV)
  plt.hlines(0, min(V), max(V), colors='black', linestyles='dashed')
  plt.vlines(-75, min(dV), max(dV), colors='black', linestyles='dashed')

  plt.text(-50, 1, 'Positive')
  plt.text(-50, -2, 'Negative')
  plt.text(E_L, max(dV) + 1, r'$E_L$')
  plt.xlabel(r'$V(t)$ (mV)')
  plt.ylabel(r'$\frac{dV}{dt}=\frac{-(V-E_L)}{\tau_m}$')
  plt.ylim(-8, 2)
  plt.show()
```

[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D4_Calculus/solutions/W0D4_Tutorial2_Solution_eab10b2d.py)

1. For $V>-75$ mV, the derivative is negative.

2. For $V<-75$ mV, the derivative is positive.

3. For $V=-75$ mV, the derivative is equal to $0$ is and a stable point when nothing changes.

### Section 2.1.1: Exact Solution of the LIF model without input
The LIF model has the exact solution:
\begin{align*}
V(t)=&\ E_L+(V_{reset}-E_L)e^{\frac{-t}{\tau_m}}\\
\end{align*}

where $\tau_m$ is the time constant, $V$ is the membrane potential,  $E_L$ is the resting potential, and $V_{reset}$ is the initial membrane potential.

 Click here for further details (from video)

Similar to the population equation, we need an initial membrane potential at time $0$ to solve the LIF model.

With this equation
\begin{align}
\frac{dV}{dt} &= \frac{-(V-E_L)}{\tau_m}\,\\
V(0)&=V_{reset},
\end{align}
where is $V_{reset}$ is called the reset potential.

The LIF model has the exact solution:
\begin{align*}
V(t)=&\ E_L+(V_{reset}-E_L)e^{\frac{-t}{\tau_m}}\\
\text{ which can be written as: }\\
\begin{matrix}\text{"Current membrane} \\ \text{potential}"\end{matrix}=&\text{"Resting potential"}+\begin{matrix}\text{"Reset potential minus resting potential} \\ \text{times exponential with rate one over time constant."}\end{matrix}\\
\end{align*}

#### Interactive Demo 2.1.1: Initial Condition $V_{reset}$
This exercise is to get an intuitive feel of how the different initial conditions $V_{reset}$ impacts the differential equation of the LIF and the exact solution for the equation:

\begin{align}
\frac{dV}{dt} &= \frac{-(V-E_L)}{\tau_m}\,\\
\end{align}
with the parameters set as:
* `E_L = -75,`
* `tau_m = 10.`

The panel on the left-hand side plot the change in membrane potential $\frac{dV}{dt}$ is a function of membrane potential $V$ and right-hand side panel plots the exact solution $V$ as a function of time $t,$ the green dot in both panels is the reset potential $V_{reset}$.

Pay close attention to when $V_{reset}=E_L=-75$mV.

1. How does the solution look with initial values of $V_{reset} < -75$?
2. How does the solution look with initial values of $V_{reset} > -75$?
3. How does the solution look with initial values of $V_{reset} = -75$?

```python
#@markdown Make sure you execute this cell to enable the widget!
my_layout.width = '450px'
@widgets.interact(
  V_reset=widgets.FloatSlider(-77., min=-91., max=-61., step=2,
                              layout=my_layout)
)

def V_reset_widget(V_reset):
  plot_V_no_input(V_reset)
```

    interactive(children=(FloatSlider(value=-77.0, description='V_reset', layout=Layout(width='450px'), max=-61.0,…

[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D4_Calculus/solutions/W0D4_Tutorial2_Solution_66938503.py)

1. Initial Values of $V_{reset} < -75$ result in the solution increasing to
  -75mV because $\frac{dV}{dt} > 0$.

2. Initial Values of $V_{reset} > -75$ result in the solution decreasing to
  -75mV because $\frac{dV}{dt}  < 0$.

3. Initial Values of $V_{reset} = -75$ result in a constant $V = -75$ mV
  because $\frac{dV}{dt}  = 0$ (Stable point).

## Section 2.2: LIF with input
We will re-introduce the input $I$ and membrane resistance $R_m$ giving the original equation:

\begin{align}
\tau_m\frac{dV}{dt} = -(V-E_L) + \color{blue}{R_mI}\,
\end{align}

The input can be other neurons or sensory information.

### Interactive Demo 2.2: The Impact of Input
The interactive plot below manipulates $I$ in the differential equation.

- With increasing input, how does the $\frac{dV}{dt}$ change? How would this impact the solution?

```python
# @markdown Make sure you execute this cell to enable the widget!
my_layout.width = '450px'
@widgets.interact(
  I=widgets.FloatSlider(3., min=0., max=20., step=2,
                                layout=my_layout)
)

def Pop_widget(I):
  plot_dVdt(I=I)
  plt.show()
```

    interactive(children=(FloatSlider(value=3.0, description='I', layout=Layout(width='450px'), max=20.0, step=2.0…

[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D4_Calculus/solutions/W0D4_Tutorial2_Solution_fba302c0.py)

dV/dt becomes bigger and less of it is below 0. This means the solution will increase well beyond what is bioligically plausible

### Section 2.2.1: LIF exact solution

The LIF with a constant input has a known exact solution:
\begin{align*}
V(t)=&\ E_L+R_mI+(V_{reset}-E_L-R_mI)e^{\frac{-t}{\tau_m}}\\
\text{which is written as:}\\
\begin{matrix}\text{"Current membrane} \\ \text{potential"}\end{matrix}=&\text{"Resting potential"}+\begin{matrix}\text{"Reset potential minus resting potential} \\ \text{times exponential with rate one over time constant." }\end{matrix}\\
\end{align*}

The plot below shows the exact solution of the membrane potential with the parameters set as:
* `V_reset = -75,`
* `E_L = -75,`
* `tau_m = 10,`
* `R_m = 10,`
* `I = 10.`

Ask yourself, does the result make biological sense? If not, what would you change? We'll delve into this in the next section

```python
# @markdown Make sure you execute this cell to see the exact solution
dt = 0.5
t_rest = 0

t = np.arange(0, 1000, dt)

tau_m = 10
R_m = 10
V_reset = E_L = -75

I = 10

V = E_L + R_m*I + (V_reset - E_L - R_m*I) * np.exp(-(t)/tau_m)

with plt.xkcd():

  fig = plt.figure(figsize=(6, 4))
  plt.plot(t,V)
  plt.ylabel('V (mV)')
  plt.xlabel('time (ms)')
  plt.show()
```

## Section 2.3: Maths is one thing, but neuroscience matters

###  Video 4: Adding firing to the LIF

```python
# @title Video 4: Adding firing to the LIF
from ipywidgets import widgets

out2 = widgets.Output()
with out2:
  from IPython.display import IFrame
  class BiliVideo(IFrame):
      def __init__(self, id, page=1, width=400, height=300, **kwargs):
          self.id=id
          src = 'https://player.bilibili.com/player.html?bvid={0}&page={1}'.format(id, page)
          super(BiliVideo, self).__init__(src, width, height, **kwargs)

  video = BiliVideo(id="BV1gX4y1P7pZ", width=730, height=410, fs=1)
  print('Video available at https://www.bilibili.com/video/{0}'.format(video.id))
  display(video)

out1 = widgets.Output()
with out1:
  from IPython.display import YouTubeVideo
  video = YouTubeVideo(id="rLQk-vXRaX0", width=730, height=410, fs=1, rel=0)
  print('Video available at https://youtube.com/watch?v=' + video.id)
  display(video)

out = widgets.Tab([out1, out2])
out.set_title(0, 'Youtube')
out.set_title(1, 'Bilibili')

display(out)
```

    Tab(children=(Output(), Output()), _titles={'0': 'Youtube', '1': 'Bilibili'})

This video first recaps the introduction of input to the leaky integrate and fire model and then delves into how we add spiking behavior (or firing) to the model.

 Click here for text recap of video

While the mathematics of the exact solution is exact, it is not biologically valid as a neuron spikes and definitely does not plateau at a very positive value.

To model the firing of a spike, we must have a threshold voltage $V_{th}$ such that if the voltage $V(t)$ goes above it, the neuron spikes
$$V(t)>V_{th}.$$
We must record the time of spike $t_{isi}$  and count the number of spikes
$$t_{isi}=t, $$
$$𝑆𝑝𝑖𝑘𝑒=𝑆𝑝𝑖𝑘𝑒+1.$$
Then reset the membrane voltage $V(t)$
$$V(t_{isi} )=V_{Reset}.$$

To take into account the spike the exact solution becomes:
\begin{align*}
V(t)=&\ E_L+R_mI+(V_{reset}-E_L-R_mI)e^{\frac{-(t-t_{isi})}{\tau_m}},&\qquad V(t)<V_{th} \\
V(t)=&V_{reset},&\qquad V(t)>V_{th}\\
Spike=&Spike+1,&\\
t_{isi}=&t,\\
\end{align*}
while this does make the neuron spike, it introduces a discontinuity which is not as elegant mathematically as it could be, but it gets results so that is good.
</detail>

### Interactive Demo 2.3.1: Input on spikes
This exercise show the relationship between firing rate and the Input for exact solution `V` of the LIF:
$$
V(t)=\ E_L+R_mI+(V_{reset}-E_L-R_mI)e^{\frac{-(t-t_{isi})}{\tau_m}},
$$
with the parameters set as:
* `V_reset = -75,`
* `E_L = -75,`
* `tau_m = 10,`
* `R_m = 10.`

Below is a figure with three panels;
* the top panel is the input, $I,$
* the middle panel is the membrane potential $V(t)$ to illustrative the spike $V(t)$ is set to $0$ and then reset to $-75$ mV,
* the bottom panel is the raster plot with each dot indicating a spike.

First, as electrophysiologist normally listen to spikes when conducting experiments, listen to the music of the firing rate for a single value of $I$. (Note the audio doesn't work in some browsers so don't worry about it if you can't hear anything)

```python
# @markdown Make sure you execute this cell to be able to hear the neuron
I = 3
t = np.arange(0, 1000, dt)
Spike, Spike_time, V = Exact_Integrate_and_Fire(I, t)

plot_IF(t, V, I, Spike_time)
ipd.Audio(V, rate=len(V))
```

Manipulate the input into the LIF to see the impact of input on the firing pattern (rate).

*  What is the effect of $I$ on spiking?
* Is this biologically valid?

```python
# @markdown Make sure you execute this cell to enable the widget!
my_layout.width = '450px'
@widgets.interact(
  I=widgets.FloatSlider(3, min=2.0, max=4., step=.1,
                              layout=my_layout)
)

def Pop_widget(I):
  Spike, Spike_time, V = Exact_Integrate_and_Fire(I, t)
  plot_IF(t, V, I, Spike_time)
```

    interactive(children=(FloatSlider(value=3.0, description='I', layout=Layout(width='450px'), max=4.0, min=2.0),…

[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D4_Calculus/solutions/W0D4_Tutorial2_Solution_f60cec6c.py)

1. As $I$ increases, the number of spikes increases.

2. No, as there is a limit to the number of spikes due to a refractory period, which is not accounted for in this model.

## Section 2.4 Firing Rate as a function of Input

The firing frequency of a neuron plotted as a function of current is called an input-output curve (F–I curve). It is also known as a transfer function, which you came across in the previous tutorial. This function is one of the starting points for the rate model, which extends from modelling single neurons to the firing rate of a collection of neurons.

By fitting this to a function, we can start to generalise the firing pattern of many neurons, which can be used to build rate models which is discussed later in Neuromatch.

```python
# @markdown *Execture this cell to visualize the FI curve*
I_range = np.arange(2.0, 4.0, 0.1)
Spike_rate = np.ones(len(I_range))

for i, I in enumerate(I_range):
  Spike_rate[i], _, _ = Exact_Integrate_and_Fire(I, t)

with plt.xkcd():
  fig = plt.figure(figsize=(6, 4))
  plt.plot(I_range,Spike_rate)
  plt.xlabel('Input Current (nA)')
  plt.ylabel('Spikes per Second (Hz)')
  plt.show()
```

The LIF model is a very nice differential equation to start with in computational neuroscience as it has been used as a building block for many papers that simulate neuronal response.

__Strengths  of LIF model:__
+ Has an exact solution;
+ Easy to interpret;
+ Great to build network of neurons.

__Weaknesses of the LIF model:__
- Spiking is a discontinuity;
- Abstraction from biology;
- Cannot generate different spiking patterns.

---
# Summary

##  Video 5: Summary

```python
# @title Video 5: Summary
from ipywidgets import widgets

out2 = widgets.Output()
with out2:
  from IPython.display import IFrame
  class BiliVideo(IFrame):
      def __init__(self, id, page=1, width=400, height=300, **kwargs):
          self.id=id
          src = 'https://player.bilibili.com/player.html?bvid={0}&page={1}'.format(id, page)
          super(BiliVideo, self).__init__(src, width, height, **kwargs)

  video = BiliVideo(id="BV1jV411x7t9", width=730, height=410, fs=1)
  print('Video available at https://www.bilibili.com/video/{0}'.format(video.id))
  display(video)

out1 = widgets.Output()
with out1:
  from IPython.display import YouTubeVideo
  video = YouTubeVideo(id="VzwLAW5p4ao", width=730, height=410, fs=1, rel=0)
  print('Video available at https://youtube.com/watch?v=' + video.id)
  display(video)

out = widgets.Tab([out1, out2])
out.set_title(0, 'Youtube')
out.set_title(1, 'Bilibili')

display(out)
```

    Tab(children=(Output(), Output()), _titles={'0': 'Youtube', '1': 'Bilibili'})

In this tutorial, we have seen two differential equations, the population differential equations and the leaky integrate and fire model.

We learned about:
* The motivation for differential equations.
* An intuitive relationship between the solution and the form of the differential equation.
* How different parameters of the differential equation impact the solution.
* The strengths and limitations of the simple differential equations.

---
# Links to Neuromatch Days

Differential equations turn up in a number of different Neuromatch days:
* The LIF model is discussed in more details in Model Types (Week 1 Day 1) and Real Neurons (Week 2 Day 3).
* Drift Diffusion model which is a differential equation for decision making is discussed in Linear Systems (Week 2 Day 2).
* Systems of differential equations are discussed in Linear Systems (Week 2 Day 2) and Dynamic Networks (Week 2 Day 4).

---
# References
1. Lotka, A. L, (1920)   Analytical note on certain rhythmic  relations inorganic  systems.Proceedings of the National Academy of Sciences,6(7):410–415,1920.

2. Brunel N, van Rossum MC. Lapicque's 1907 paper: from frogs to integrate-and-fire. Biol Cybern. 2007 Dec;97(5-6):337-9. doi: 10.1007/s00422-007-0190-0. Epub 2007 Oct 30. PMID: 17968583.

# Bibliography
1. Dayan, P., & Abbott, L. F. (2001). Theoretical neuroscience: computational and mathematical modeling of neural systems. Computational Neuroscience Series.
2. Strogatz, S. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity), Westview Press; 2 edition (29 July 2014)

## Supplemental Popular Reading List
1. Lindsay, G. (2021). Models of the Mind: How Physics, Engineering and Mathematics Have Shaped Our Understanding of the Brain. Bloomsbury Publishing.
2. Strogatz, S. (2004). Sync: The emerging science of spontaneous order. Penguin UK.

## Popular Podcast
1. Strogatz, S. (Host). (2020-), Joy of X  https://www.quantamagazine.org/tag/the-joy-of-x/ Quanta Magazine

ext: ipynb


text:
# Домашнее задание "Математический анализ. Производная".

### Уровень 0:

Посчитайте значение первой и второй производной функции:

$f(x) = x^5 + 4sin(2x) + cos(3x+3)$

В точке x = 1.

Analytical solution:

$f'(x) = 5x^4 + 8cos(2x) - 3sin(3x+3)$

$f''(x) = 20x^3 - 16sin(2x) - 9cos(3x+3)$

```python
import math
from scipy.misc import derivative

def f(x):
    y = x**5+4*math.sin(2*x)+math.cos(3*x+3)
    return y

def f1(x):
    y = 5*x**4+8*math.cos(2*x)-3*math.sin(3*x+3)
    return y

def f2(x):
    y = 20*x**3-16*math.sin(2*x)-9*math.cos(3*x+3)
    return y

dx1 = derivative(f, 1.0, dx=1e-6)
print("dx1: {}; f1(1):{}".format(dx1, f1(1.0)))
dx2 = derivative(f, 1.0, dx=1e-6, n=2)
print("dx2: {}; f2(1):{}".format(dx2, f2(1.0)))
```

    dx1: 2.509071801881646; f1(1):2.5090718022196383
    dx2: -3.19033688356285; f2(1):-3.1902914090642014

### Уровень 1:

Постройте график первой и второй производной функции:

$f(x) = sin^5(2x+1)$

На отрезке [-5; 5].

```python
from sympy import *
from numpy import linspace
import matplotlib.pyplot as pl

x = symbols('x')
y = sin(2 * x + 1) ** 5
```

```python
y1 = diff(y, x)
print(format(y1))

lam_y1 = lambdify(x, y1, modules=['numpy'])

x_vals = linspace(-5, 5, 1000)
y_vals = lam_y1(x_vals)

pl.plot(x_vals, y_vals)
pl.show()
```

```python
y2 = diff(y1, x)
print(format(y2))

lam_y2 = lambdify(x, y2, modules=['numpy'])

x_vals = linspace(-5, 5, 1000)
y_vals = lam_y2(x_vals)

pl.plot(x_vals, y_vals)
pl.show()
```

```python

```