<jupyter_start><jupyter_text>[![Kaggle](https://kaggle.com/static/images/open-in-kaggle.svg)](https://kaggle.com/kernels/welcome?src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W0D3_LinearAlgebra/W0D3_Tutorial3.ipynb)  Bonus Tutorial: Discrete Dynamical Systems**Week 0, Day 3: Linear Algebra****By Neuromatch Academy**__Content creators:__ <NAME>__Production editors:__ <NAME> **Our 2021 Sponsors, including Presenting Sponsor Facebook Reality Labs** --- Tutorial ObjectivesIn this tutorial, we will start to gain an intuition for how eigenvalues and eigenvectors can be helpful for understanding dynamical systems. We will focus on a discrete dynamical system consisting of two neurons. By the end of the tutorial, you will:* Predict whether the firing rates of interconnected model neurons will explode or decay based on the eigenvalues of the weight matrix.* Apply ideas from previous tutorials (linear combination, basis vectors, etc) to understand a new concept --- Setup<jupyter_code># Imports

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
#@title Figure settings
import ipywidgets as widgets       # interactive display
%config InlineBackend.figure_format = 'retina'
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")
#@title Plotting functions

def plot_circuit_responses(u, W, eigenstuff = False, xlim='default', ylim='default'):
    fig, ax = plt.subplots(1, 1, figsize=(10,10))

    # Set up axis limits
    if xlim =='default':
      extreme = np.maximum(np.abs(np.min(u)), np.max(u))
      xlim = [- extreme, extreme]
    if ylim == 'default':
      extreme = np.maximum(np.abs(np.min(u)), np.max(u))
      ylim = [- extreme, extreme]

    # Set up look
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    cs = plt.rcParams['axes.prop_cycle'].by_key()['color']*10
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)

    # Set up tracking textz
    tracker_text = ax.text(.5, .9, "", color='w', fontsize=20, verticalalignment='top', horizontalalignment='left', transform=ax.transAxes)

    # Plot eigenvectors
    if eigenstuff:
      eigvals, eigvecs = np.linalg.eig(W)

      if np.abs(eigvals[0]) < np.abs(eigvals[1]):
        lc1 = 'c'
        lc2 = 'g'
      else:
        lc1 = 'g'
        lc2 = 'c'

      ax.plot(np.arange(-10000, 10000)*eigvecs[0, 0], np.arange(-10000, 10000)*eigvecs[1, 0],lc1, alpha=.5, label = r'$\mathbf{v}_1$')
      ax.plot(np.arange(-10000, 10000)*eigvecs[0, 1], np.arange(-10000, 10000)*eigvecs[1, 1], lc2, alpha=.5, label = r'$\mathbf{v}_2$')

      ax.legend()

    # Set up scatter
    cmap = plt.cm.coolwarm
    norm = plt.Normalize(vmin=0, vmax=u.shape[1])
    ax.plot(u[0,:], u[1, :], alpha=.4, zorder = 0)
    scatter = ax.scatter(u[0, :], u[1, :], alpha=1, c = cmap(norm(np.arange(u.shape[1]))), zorder = 1)


    ax.set(xlabel = 'Neuron 1 Firing Rate ($a_1$)', ylabel = 'Neuron 2 Firing Rate ($a_2$)', title = 'Neural firing over time')

    fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap),
                ax=ax, label = 'Time step')
#@title Helper functions

def get_eigval_specified_matrix(target_eig):
  """Generates matrix with specified eigvals

    Args:
      target_eig (list): list of target eigenvalues, can be real or complex,
              should be length 2 unless you desire repeated eigenvalues
              with the same eigenvector, in which case length 1

    Returns:
      ndarray: 2 x 2 matrix with target eigvals

  """

  # Set up two eigenvectors
  V = np.array([[1, 1], [-1, 1]]).astype('float')
  for i in range(2):
    V[:,i] = V[:,i]/np.linalg.norm(V[:,i])

  # Get matrix with target eigenvalues
  if type(target_eig[0]) == int or type(target_eig[0]) == float:

    if len(target_eig) == 2: # distinct eigvecs (not necessarily distinct eigvals)

        D = np.diag(target_eig)
        A = V @ D @ np.linalg.inv(V)

    else: # repeated with same vec
      summed = 2*target_eig[0]

      a = summed-3
      d = 3
      bc = target_eig[0]**2 - a*d
      factors = [n for n in range(1, bc+ 1) if bc % n == 0]
      b = factors[int(np.floor(len(factors)/2))]
      c = bc/-b

      A = np.array([[a, b], [c, d]])

  elif type(target_eig[0]) == complex:

      C = [np.real(V[:,0]), np.real(V[:,1])]
      B = np.array([[np.real(target_eig[0]), np.imag(target_eig[0])], [-np.imag(target_eig[0]), np.real(target_eig[0])]]).squeeze()
      A = C @ B @ np.linalg.inv(C)

  return A<jupyter_output><empty_output><jupyter_text>--- Section 1: Defining a neural circuit **Please note a correction to the video appearing below it**<jupyter_code># @title Video 1: A neural circuit
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="BV1564y1b78X", width=854, height=480, 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="XJu2bQgnn9M", width=854, height=480, 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)<jupyter_output><empty_output><jupyter_text>This video covers our neural circuit we will be working with in this tutorial.**Please note this correction to the video: I had wrong subscripts for weights in the equations and matrix. These should be:**$$a_{1, t} = w_{1, 1}a_{1, t-1} + w_{2, 1}a_{2, t-1} $$$$a_{2, t} = w_{1, 2}a_{1, t-1} + w_{2, 2}a_{2, t-1} $$$$\mathbf{W} = \begin{bmatrix} w_{1, 1} & w_{2, 1} \\ w_{1, 2} & w_{2, 2} \end{bmatrix}, \mathbf{a}_{t} =   \begin{bmatrix} a_{1, t} \\ a_{2, t} \end{bmatrix}$$ Click here for text recap of video In previous tutorials, we have looked at static models of LGN neurons based on the responses of retinal neurons. By static, I mean that we just looked at a single time point.Let's now introduce the concept of time. We will chop time up into little bins and look at the activity of neurons in each bin. That is, we will work in a **discrete** time framework. For example, if each bin is 1 second long, we will look at the firing rate of each neuron at intervals of a second.Instead of examining pre- and post- synaptic neurons, we will examine two neurons in one area that are connected. In our model, the activity of neuron 1 at one time bin depends on the activities of both neurons during the previous time bin multiplied by the respective weights from itself and neuron 2. It might seem weird for a neuron to have a weight to itself - this is abstracting away some biological details but basically conveys how much the neural activity depends on its history. (Throughout this course, we'll see lots of neuron models and how some model biological detail more faithfully while others abstract.)We will refer to the activity of neuron i during time bin j as $a_{i, j}$. The weight from neuron y to neuron x will be $w_{y, x}$. With this helpful notation, we can write an equation for the activity of neuron 1 at time bin t:$$a_{1, t} = w_{1, 1}a_{1, t-1} + w_{2, 1}a_{2, t-1} $$And the symmetric model is true of neuron 2:$$a_{2, t} = w_{1, 2}a_{1, t-1} + w_{2, 2}a_{2, t-1} $$This is already a mess of subscript numbers - luckily we can use matrices and vectors once again and our model becomes: $$\mathbf{a}_{t} = \mathbf{W}\mathbf{a}_{t-1} $$where:$$\mathbf{W} = \begin{bmatrix} w_{1, 1} & w_{2, 1} \\ w_{1, 2} & w_{2, 2} \end{bmatrix}, \mathbf{a}_{t} =   \begin{bmatrix} a_{1, t} \\ a_{2, t} \end{bmatrix}$$It turns out that this is a **discrete dynamical system**. Dynamical systems are concerned with how quantities evolve other time, in this case our neural firing rates. When we model the evolution of quantities over time using a discrete time framework, it is, unsurprisingly, a discrete dynamical system. We will see continuous dynamical systems (where we embrace the full continuity of time) tomorrow and later in the comp neuro course during W2D2: Linear Dynamics.  Coding Exercise 1: Implementing the circuitIn this exercise, you will implement the function `circuit_implementation`. Given a weight matrix, initial activities at time 0, and a number of time bins to model, this function calculates the neural firing rates at each time bin.We will use initial firing rates of 1 for both neurons:$$\mathbf{a}_0 = \begin{bmatrix}1 \\1  \\\end{bmatrix}$$and the weight matrix:$$\mathbf{W} = \begin{bmatrix} 1 & 0.2 \\0.1 & 1 \\ \end{bmatrix}$$We will look at activity over 30 time steps. As before, we will allow our firing rates to be negative, despite this not being possible biologically.<jupyter_code>def circuit_implementation(W, a0, T):
  """ Simulate the responses of N neurons over time given their connections

  Args:
    W (ndarray): weight matrix of synaptic connections, should be N x N
    a0 (ndarray): initial condition or input vector, should be N,
    T (scalar): number of time steps to run simulation for

  Returns:
    a (ndarray): the neural responses over time, should be N x T

  """

  # Compute the number of neurons
  N = W.shape[0]

  # Initialize empty response array and initial condition
  a = np.zeros((N, T))
  a[:, 0]  = a0

  #################################################
  ## TODO for students ##
  # Fill out function and remove
  raise NotImplementedError("Student exercise: Complete circuit_implementation")
  #################################################

  # Loop over time steps and compute u(t+1)
  for i_t in range(1, T):
      a[:, i_t] = ...

  return a


# Define W, u0, T
W = np.array([[1, .2], [.1, 1]])
a0 = np.array([1, 1])
T = 30

# Get neural activities
a = circuit_implementation(W, a0, T)

# Visualize neural activities
plot_circuit_responses(a, W)
# to_remove solution
def circuit_implementation(W, a0, T):
  """ Simulate the responses of N neurons over time given their connections

  Args:
    W (ndarray): weight matrix of synaptic connections, should be N x N
    a0 (ndarray): initial condition or input vector, should be N,
    T (scalar): number of time steps to run simulation for

  Returns:
    a (ndarray): the neural responses over time, should be N x T

  """

  # Compute the number of neurons
  N = W.shape[0]

  # Initialize empty response array and initial condition
  a = np.zeros((N, T))
  a[:, 0]  = a0

  # Loop over time steps and compute u(t+1)
  for i_t in range(1, T):
      a[:, i_t] = W @ a[:, i_t-1]

  return a


# Define W, u0, T
W = np.array([[1, .2], [.1, 1]])
a0 = np.array([1, 1])
T = 30

# Get neural activities
a = circuit_implementation(W, a0, T)

# Visualize neural activities
with plt.xkcd():
  plot_circuit_responses(a, W)<jupyter_output><empty_output><jupyter_text>The firing rates of both neurons are exploding to infinity over time. Let's now see what happens with a different weight matrix:$$\mathbf{W} = \begin{bmatrix} 0.2 & 0.1 \\1 & 0.2 \\ \end{bmatrix}$$<jupyter_code># @markdown Execute this cell to visualize activity over time

# Define W, a0, T
W = np.array([[.2, .1], [1, .2]])
a0 = np.array([1, 1])
T = 30

# Get neural activities
a = circuit_implementation(W, a0, T)

# Visualize neural activities
plot_circuit_responses(a, W)<jupyter_output><empty_output><jupyter_text>We can see that with this weight matrix, the firing rates are decaying towards zero. It turns out that we could have predicted this by looking at the eigenvalues of the weight matrices, as we'll see in the next section. --- Section 2: Understanding dynamics using eigenstuffAs we'll see in this section, eigenvectors and eigenvalues are incredibly useful for understanding the evolution of the neural firing rates, and discrete dynamical systems in general.  Section 2.1: Initial firing rates along an eigenvector **Please note this correction to the video: the equation at the bottom of slide 115 (at 2:31) should be: $a_i = W^ia_0$**<jupyter_code># @title Video 2: Looking at activity along an eigenvector
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="BV1Tf4y1b7de", width=854, height=480, 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="BnSvcWbYf8g", width=854, height=480, 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)<jupyter_output><empty_output><jupyter_text>This video covers what happens if our initial firing rates fall along one of the eigenvectors of the weight matrix. Click here for text recap of video **Rewriting our circuit equation**In our neural circuit, we are modeling the activities at a time step as:$$\mathbf{a}_{t} = \mathbf{W}\mathbf{a}_{t-1} $$Let's start at time  step 1:$$\mathbf{a}_{1} = \mathbf{W}\mathbf{a}_{0} $$And move on to time step 2:$$\mathbf{a}_{2} = \mathbf{W}\mathbf{a}_{1} $$In the above equation, we can subsitute in $\mathbf{a}_{1} = \mathbf{W}\mathbf{a}_{0}$:$$\mathbf{a}_{2} = \mathbf{W}\mathbf{W}\mathbf{a}_{0} = \mathbf{W}^2 \mathbf{a}_{0}$$We can keep doing this with subsequent time steps:$$\mathbf{a}_{3} = \mathbf{W}\mathbf{a}_{2} = \mathbf{W}\mathbf{W}^2 \mathbf{a}_{0} = \mathbf{W}^3\mathbf{a}_{0}  $$$$\mathbf{a}_{4} = \mathbf{W}\mathbf{a}_{3} = \mathbf{W}\mathbf{W}^3 \mathbf{a}_{0} = \mathbf{W}^4\mathbf{a}_{0}  $$This means that we can write the activity at any point as:$$\mathbf{a}_{i} = \mathbf{W}^i\mathbf{a}_{0}  $$**Initial firing rates along an eigenvector**Remember from the last tutorial, that an eigenvector of matrix $\mathbf{W}$ is a vector that becomes a scalar multiple (eigenvalue) of itself when multiplied by that matrix:$$\mathbf{W}\mathbf{v} = \lambda\mathbf{v}$$Let's look at what happens if the initial firing rates in our neural circuit lie along that eigenvector, using the same substitution method as in the previous section:$$\mathbf{a}_{0} = \mathbf{v} $$$$\mathbf{a}_{1} = \mathbf{W}\mathbf{a}_0 = \mathbf{W}\mathbf{v} = \lambda\mathbf{v} $$$$\mathbf{a}_{2} = \mathbf{W}\mathbf{a}_1 = \mathbf{W}\lambda\mathbf{v} = \lambda\mathbf{W}\mathbf{v} = \lambda^2\mathbf{v}$$$$\mathbf{a}_{3} = \mathbf{W}\mathbf{a}_2 = \mathbf{W}\lambda^2\mathbf{v} = \lambda^2\mathbf{W}\mathbf{v} = \lambda^3\mathbf{v}$$$$...$$$$\mathbf{a}_i = \lambda^i\mathbf{v}$$The activities at any time step equal a scalar times the initial activities. In other words, if the initial activities lie along an eigenvector, the activities will only evolve along that eigenvector.   Interactive demo 2.1: Changing the eigenvalueLet's visualize what happens if the initial activities of the neurons lie along an eigenvector and think about how this depends on the eigenvalue.The interactive demo below is the same visualization you saw in Section 1, but now we also plot the eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$.Questions:1.  What happens if the eigenvalue is large (2)?2.  What happens if you move the eigenvalue from 2 to towards 0? 3.  What happens with negative eigenvalues?<jupyter_code># @markdown Execute this cell to enable the widget. Please be patient for a few seconds after you change the slider

@widgets.interact(eigenvalue = widgets.FloatSlider(value=0.5, min=-2, max=2, step=0.2))
def plot_system(eigenvalue):

  # Get weight matrix with specified eigenvalues
  W = get_eigval_specified_matrix([eigenvalue, eigenvalue])

  # Get initial condition
  u0 = np.array([1, 1])

  # Get neural activities
  u = circuit_implementation(W, u0, 10)

  # Visualize neural activities
  plot_circuit_responses(u, W, eigenstuff = True, xlim = [-15, 15], ylim = [-15, 15])
# to_remove explanation

# 1) With the eigenvalue = 2, the activities of the neurons explode towards infinity, along
#.   the eigenvector.

# 2) At eigenvalue = 1, there is a shift in what happens. With the eigenvalue above 1,
#.  the activites always explode. Once the eigenvalue is below 1, the activities decay to 0.
#.  If the eigenvalue equals 1, the activities never differ from the initial condition.
#.  This makes sense with the equation above. Lambda is raised to a power when computing activities:
#.  if it's a fraction, this term will get smaller so the activities will. If above 1, this term
#.   will explore so the activities will.

# 3) If the eigenvalue is between -1 and 0, the neural activities jump across the
#.   origin repeatedly along the eigenvector but eventually decay to 0. If the eigenvalue is below -1, the
#.   activities jump across the origin repeatedly along the eigenvector but explode to
#.   positive or negative infinity. Once again, this makes sense if you think through the equation above.<jupyter_output><empty_output><jupyter_text>Section 2.2: Understanding general dynamics using eigenstuff<jupyter_code># @title Video 3: Understanding general dynamics using eigenstuff
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="BV1fX4y1c7y7", width=854, height=480, 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="66DExLQPzPI", width=854, height=480, 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)<jupyter_output><empty_output><jupyter_text>This video covers how the eigenvalues and eigenvectors of the weight matrix affect the neural activities more generally. Click here for text recap of video We now know that if our initial activities (or initial condition) fall on an eigenvector of $\mathbf{W}$, the activities will evolve along that line, either exploding to infinity if the absolute value of the eigenvalue is above 1 or decaying to the origin it it is below 1. What if our initial condition doesn't fall along the eigenvector though?To understand what will happen, we will use the ideas of basis vectors and linear combinations from Tutorial 1.Let's assume for now that our weight matrix has two distinct eigenvectors ($\mathbf{v}_1$ and $\mathbf{v}_2$) with corresponding eigenvalues $\lambda_1$ and $\lambda_2$, and that these eigenvectors form a basis for 2D space. That means we can write any vector in 2D space as a linear combination of our eigenvectors, including our initial activity vector:$$\mathbf{a}_0 = c_1\mathbf{v}_1 + c_2\mathbf{v}_2 $$Let's compute the next time step, using our previous strategy of substitution:$$\begin{align}\mathbf{a}_1 &= \mathbf{W}\mathbf{a}_0\\ &= \mathbf{W}(c_1\mathbf{v}_1 + c_2\mathbf{v}_2) \\ &= c_1\mathbf{W}\mathbf{v}_1 + c_2\mathbf{W}\mathbf{v}_2 \\ &= c_1\lambda_1\mathbf{v}_1 + c_2\lambda_2\mathbf{v}_2 \end{align} $$All activities can be written as:$$\mathbf{a}_i = c_1\lambda_1^i\mathbf{v}_1 + c_2\lambda_2^i\mathbf{v}_2 $$ We'll see what this means for our system in the next demo. Before we get to the neural circuit, refresh your memory on linear combinations briefly by looking at our widget from tutorial 1 below. What happens when the absolute values of both scalar multiples of both vectors is big? What about when one is big and one is small?<jupyter_code>#@markdown Make sure you execute this cell to enable the widget! Move the sliders for “a” and “b”. After releasing the slider, be patient for a couple of seconds to see the desired change.


def plot_arrows(x, y, a_times_x, b_times_y):
  fig, ax = plt.subplots(figsize=(10, 7))

  ax.spines['top'].set_color('none')
  ax.spines['bottom'].set_position('zero')
  ax.spines['left'].set_position('zero')
  ax.spines['right'].set_color('none')

  ax.set_aspect('equal', adjustable='box')
  ax.set(xlim = [-10, 10], ylim = [-10, 10], xticks = np.arange(-10, 10), yticks = np.arange(-10, 10), xticklabels = [], yticklabels = [])

  ax.grid(alpha=.4)

  z = a_times_x + b_times_y
  z_arr = ax.arrow(0, 0, z[0], z[1], width=.08, color='k', length_includes_head = True);

  x_orig, = ax.plot([0, x[0]], [0, x[1]], '--', color='#648FFF')
  y_orig, = ax.plot([0, y[0]], [0, y[1]], '--',  color='#DC267F')

  ax_arr = ax.arrow(0, 0, a_times_x[0], a_times_x[1], width=.08, color='#648FFF', length_includes_head = True);
  by_arr = ax.arrow(0, 0, b_times_y[0], b_times_y[1], width=.08, color='#DC267F', length_includes_head = True);

  ax.plot([a_times_x[0], z[0]], [a_times_x[1], z[1]], '--k')
  ax.plot([b_times_y[0], z[0]], [b_times_y[1], z[1]], '--k')


  leg = ax.legend([x_orig, y_orig, ax_arr, by_arr, z_arr], [r"$\mathbf{x}$", r"$\mathbf{y}$", r"$a\mathbf{x}$", r"$b\mathbf{y}$", r"$\mathbf{z} = a\mathbf{x} + b\mathbf{y}$"], handlelength = 2, fontsize = 25, loc = 'center left', bbox_to_anchor=(1.05, .5))
  for handle, label in zip(leg.legendHandles, leg.texts):
      try:
          label.set_color(handle.get_facecolor())
      except:
          label.set_color(handle.get_color())
      #handle.set_visible(False)

@widgets.interact(a = widgets.FloatSlider(value=1.0, min=-4, max=4, step=0.1), b = widgets.FloatSlider(value=1.0, min=-4, max=4, step=0.1))
def plot_linear_combination(a, b):
  x = np.array([3, 1])
  y = np.array([-1, 2])

  plot_arrows(x, y, a*x, b*y)<jupyter_output><empty_output><jupyter_text>Interactive demo 2.2: Changing both eigenvaluesIn the demo below, you can now change both eigenvalues and the initial condition (with `a0_1` setting neuron 1 initial activity and `a0_2` setting neuron 2 initial activity). We will only look at positive eigenvalues to keep things a little more simple. We also make sure the second eigenvalue is always the smaller one (just for better visualization purposes).Think each of the following questions through based on the equation we just arrived at and then play with the demo to see if you are correct.$$\mathbf{a}_i = c_1\lambda_1^i\mathbf{v}_1 + c_2\lambda_2^i\mathbf{v}_2 $$1.  What will happen when both eigenvalues are greater than 1? Does this depend on initial condition?  Set eigenvalue1 to 2 and eigenvalue2 to 1.2 and try out different initial conditions. What do you see? 2.  What will happen when both eigenvalues are less than 1?3. What happens if one eigenvalue is below 1 and the other is above 1?<jupyter_code># @markdown Execute this cell to enable the widget (there is a small lag so be patient after changing sliders)
eigenvalue1 = widgets.FloatSlider(value=0.5, min=0.2, max=2, step=0.2)
eigenvalue2 = widgets.FloatSlider(value=0.2, min=0.2, max=0.5, step=0.2)
a0_1 = widgets.FloatSlider(value=1, min=-5, max=5, step=0.2)
a0_2 = widgets.FloatSlider(value=2, min=-5, max=5, step=0.2)

def update_range(*args):
    eigenvalue2.max = eigenvalue1.value - 0.2
eigenvalue1.observe(update_range, 'value')


def plot_system(eigenvalue1, eigenvalue2, a0_1, a0_2):

  # Get initial condition
  a0 = np.array([a0_1, a0_2])

  # Get weight matrix with specified eigenvalues
  W = get_eigval_specified_matrix([eigenvalue1, eigenvalue2])

  # Get neural activities
  u = circuit_implementation(W, a0, 10)

  # Visualize neural activities
  plot_circuit_responses(u, W, eigenstuff = True, xlim = [-15, 15], ylim = [-15, 15])

widgets.interact(plot_system, eigenvalue1 = eigenvalue1, eigenvalue2 = eigenvalue2, a0_1 = a0_1, a0_2 = a0_2)
# to_remove explanation

# 1) If both eigenvalues are above 1, the neural activities will eventually explode
#.   to infinity or negative infinity, depending on initial conditions. Tthe exact trajectory
#.   is drawn towards the eigenvector with the larger eigenvalue. This is because the larger eigenvalue
#.   will increasingly dominate the other one as it is raised to increasingly larger powers.

# 2) If both eigenvalues are below 1, the neural activity will eventually decay to 0.

#. 3) The activities will eventually explode to positive or negative infinity, unless
#.   the initial condition lies exactly on the eigenvector with the small eigenvalue. If the
#.   initial condition is near to that eigenvector, the trajectory will first go towards
#.   the origin before exploding.<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12926732673267327, 'fraction_numerical': 0.02803960396039604, 'mean_word_length': 4.3441269841269845, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 39, 'https://': 13, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '89308', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 8664, 'n_tokens_neox': 8014, 'n_words': 2911, 'path': 'tutorials/W0D3_LinearAlgebra/W0D3_Tutorial3.ipynb', 'repo_name': 'ofou/course-content', 'size': '25258'}