added EI notebook

README.md

@@ -3,7 +3,7 @@
 ## Introduction
 This package provides an implementation of a recurrent neural network simulator with NUMBA.
 The network can have multiple neural populations, different connectivity profiles (all to all, sparse, tuned, ...).
-For more info look at the config files in ./conf
+For more info look at the config files in ./conf/.
 
 ## Installation
 Provide clear instructions on how to get your development environment running.
@@ -25,7 +25,9 @@ model = Network(config_file_name, output_file_name, path_to_repo, **kwargs)
 # run the model
 model.run()
 ```
-
+There are two configs here:
+- The first one is config_bump.py which is a continuous 1 population bump attractor model as in the NB stim paper.
+- The second is config_EI.py which are standard parameters for a tuned bump attractor balance network with 2 populations.
 
 ## Contributing
 Feel free to contribute.

conf/config_EI.yml

@@ -105,7 +105,7 @@ M0: 1.0
 
 # To add an attentional switch
 # if BUMP_SWITCH[i] == 1 it sets Iext[i] to zero before stimulus presentation
-BUMP_SWITCH: [0, 0]
+BUMP_SWITCH: [1, 1]
 
 ########################
 # Plasticity

org/EI_bal.org

@@ -0,0 +1,213 @@
+#+STARTUP: fold
+#+TITLE: EI Bump Attractor Network Model
+#+PROPERTY: header-args:ipython :results both :exports both :async yes :session dual_data :kernel dual_data
+
+* notebook settings
+
+#+begin_src ipython
+  %load_ext autoreload
+  %autoreload 2
+  %reload_ext autoreload
+
+  %run ../notebooks/setup.py
+  %matplotlib inline
+  %config InlineBackend.figure_format = 'png'
+#+end_src
+
+#+RESULTS:
+: The autoreload extension is already loaded. To reload it, use:
+:   %reload_ext autoreload
+: Python exe
+: /home/leon/mambaforge/envs/dual_data/bin/python
+
+* EI Balanced Attractor Model
+Here I implemented a balanced EI network that exhibits attractor dynamics.
+** imports
+#+begin_src ipython
+  import sys
+  sys.path.insert(0, '/home/leon/tmp/rnn_numba') # put here the path to the repo
+  from src.model.rate_model import Network
+#+end_src
+
+#+RESULTS:
+
+** Single trial
+*** Simulation
+To run a simulation, first we need to define a network model.
+The class Network takes three mandatory arguments:
+
+                1. The name of the configuration file that defines the model,
+                   eg 'config_EI.yml', these files are in ../conf/ and well detailed.
+
+                2. The name of the output file that will contain the simulation data.
+                   eg 'bump'. Simulation results will be saved as a data frame to '../data/simul/bump.h5'.
+
+                3. The path to the root of this repository.
+
+One can also pass extra arguments to Network, basically any parameter that is in the config file so that it will be overwritten.
+
+#+begin_src ipython
+  REPO_ROOT = "/home/leon/tmp/rnn_numba"
+  model = Network('config_EI.yml', 'bump', REPO_ROOT, VERBOSE=1)
+#+end_src
+
+#+RESULTS:
+: Loading config from /home/leon/tmp/rnn_numba/conf/config_EI.yml
+: Saving data to /home/leon/tmp/rnn_numba/data/simul/bump.h5
+: Jab [[ 1.  -1.5]
+:  [ 1.  -1. ]]
+: Tuning, KAPPA [5. 0. 0. 0.]
+: Asymmetry, SIGMA [0. 0. 0. 0.]
+: Iext [0.5  0.25]
+
+Then one just runs the model with
+#+begin_src ipython
+  model.run()
+#+end_src
+
+#+RESULTS:
+#+begin_example
+  Generating matrix Cij
+  sparse connectivity
+  with spec cosine structure
+  sparse connectivity
+  sparse connectivity
+  sparse connectivity
+  Saving matrix to /home/leon/tmp/rnn_numba/data/matrix/Cij.npy
+  Parameters:
+  N 10000 Na [7500 2500] K 500.0 Ka [500. 500.]
+  Iext [11.18033989  5.59016994] Jab [ 0.04472136 -0.06708204  0.04472136 -0.04472136]
+  Tuning, KAPPA [5. 0. 0. 0.]
+  Asymmetry, SIGMA [0. 0. 0. 0.]
+  MF Rates: [0.25 0.5 ]
+  Transfert Func Sigmoid
+  Running simulation
+  times (s) 0.25 rates (Hz) [0.03, 0.12]
+  times (s) 0.5 rates (Hz) [0.03, 0.12]
+  times (s) 0.75 rates (Hz) [0.03, 0.11]
+  times (s) 1.0 rates (Hz) [0.03, 0.12]
+  STIM ON
+  times (s) 1.25 rates (Hz) [0.53, 0.74]
+  times (s) 1.5 rates (Hz) [0.53, 0.74]
+  STIM OFF
+  times (s) 1.75 rates (Hz) [0.36, 0.6]
+  times (s) 2.0 rates (Hz) [0.36, 0.6]
+  times (s) 2.25 rates (Hz) [0.36, 0.6]
+  times (s) 2.5 rates (Hz) [0.36, 0.6]
+  times (s) 2.75 rates (Hz) [0.36, 0.6]
+  times (s) 3.0 rates (Hz) [0.36, 0.6]
+  times (s) 3.25 rates (Hz) [0.36, 0.6]
+  times (s) 3.5 rates (Hz) [0.36, 0.6]
+  times (s) 3.75 rates (Hz) [0.36, 0.6]
+  times (s) 4.0 rates (Hz) [0.36, 0.6]
+  saving data to /home/leon/tmp/rnn_numba/data/simul/bump.h5
+  Elapsed (with compilation) = 59.40574227599427s
+#+end_example
+
+*** Analysis
+**** Imports
+#+begin_src ipython
+  import pandas as pd
+  from src.analysis.decode import decode_bump
+#+end_src
+
+#+RESULTS:
+
+**** Load data
+#+begin_src ipython
+  df = pd.read_hdf(REPO_ROOT + "/data/simul/bump.h5", mode="r")
+  df_E = df[df.neurons<7500]
+  df_I = df[df.neurons>=7500]
+
+  print(df.head())
+#+end_src
+
+#+RESULTS:
+:       rates        ff       h_E       h_I  neurons   time
+: 0  0.016484  0.516858  0.660207 -4.258005        0  0.249
+: 1  0.010508  0.068475  0.606468 -4.101016        1  0.249
+: 2  0.039072  0.148301  0.625659 -3.540477        2  0.249
+: 3  0.033662  1.185190  0.600276 -3.787884        3  0.249
+: 4  0.044454 -0.488017  0.634493 -3.749277        4  0.249
+**** Rates
+***** raster
+#+begin_src ipython
+  fig, ax = plt.subplots(1, 2, figsize=[2*width, height])
+  pt = pd.pivot_table(df, values="rates", index=["neurons"], columns="time")
+
+  sns.heatmap(pt[:7500], cmap="jet", ax=ax[0], vmax=1, vmin=0)
+  ax[0].set_yticks([0, 2500, 5000, 7500], [0, 2500, 5000, 7500])
+  ax[0].set_xticks([0, 2, 4, 6, 8], [0, 1, 2, 3, 4])
+  ax[0].set_title('Excitatory')
+
+  sns.heatmap(pt[7500:], cmap="jet", ax=ax[1], vmax=1, vmin=0)
+  ax[1].set_yticks([0, 625, 1250, 1875, 2500], [0, 625, 1250, 1875, 2500])
+  ax[1].set_xticks([0, 2, 4, 6, 8], [0, 1, 2, 3, 4])
+  ax[1].set_title('Inhibitory')
+
+  plt.show()
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/f33e2bf59d4004bf01d9b42395a2fe97b381e362.png]]
+
+***** histograms
+
+#+begin_src ipython
+  mean_df_E = df_E.groupby("neurons").mean()
+  mean_df_E[mean_df_E.rates<.001] = np.nan
+
+  mean_df_I = df_I.groupby("neurons").mean()
+  mean_df_I[mean_df_I.rates<.001] = np.nan
+
+  sns.histplot(mean_df_E, x=mean_df_E.rates, kde=True, color='r', label='E')
+  sns.histplot(mean_df_I, x=mean_df_I.rates, kde=True, color='b', label='I')
+  plt.legend(fontsize=12)
+  plt.xlabel("Rates (au)")
+  plt.show()
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/66e9897b32db80297981892d9c7c48cbb08188f0.png]]
+
+**** Tuning
+
+#+begin_src ipython
+  data = df_E.groupby(['time'])['rates'].apply(decode_bump).reset_index()
+  data[['m0', 'm1', 'phase']] = pd.DataFrame(data['rates'].tolist(), index=data.index)
+  data = data.drop(columns=['rates'])
+
+  print(data.head())
+#+end_src
+
+#+RESULTS:
+:     time        m0        m1     phase
+: 0  0.249  0.026883  0.000102  6.019133
+: 1  0.499  0.026896  0.000113  3.882348
+: 2  0.749  0.026606  0.000504  0.829742
+: 3  0.999  0.027251  0.000305  2.127723
+: 4  1.249  0.534482  0.597531  3.139320
+
+#+begin_src ipython
+  fig, ax = plt.subplots(1, 3, figsize=[2*width, height])
+
+  sns.lineplot(data=data, x='time', y='m0', legend=False, lw=2, ax=ax[0])
+  ax[0].set_xlabel('Time (s)')
+  ax[0].set_ylabel('$\mathcal{F}_0 (Hz)$')
+  ax[1].set_xticks([0, 1, 2, 3, 4])
+
+  sns.lineplot(x=data['time'], y=data['m1']/data['m0'], legend=False, lw=2, ax=ax[1])
+  ax[1].set_xlabel('Time (s)')
+  ax[1].set_ylabel('$\mathcal{F}_1 / \mathcal{F}_0$')
+  ax[1].set_xticks([0, 1, 2, 3, 4])
+
+  sns.lineplot(x=data['time'], y=data['phase']*180/np.pi, legend=False, lw=2, ax=ax[2])
+  ax[2].set_xlabel('Time (s)')
+  ax[2].set_ylabel('$\phi$ (°)')
+  ax[2].set_xticks([0, 1, 2, 3, 4])
+  ax[2].set_yticks([0, 90, 180, 270, 360])
+  plt.show()
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/2596d8031fe7e25fbeb082a671606fd6f2681fdc.png]]

org/bump.org

@@ -21,6 +21,9 @@
 : /home/leon/mambaforge/envs/dual_data/bin/python
 
 * Continuous Bump Attractor Model
+
+Here I will explain how to use the package to run a one population continuous attractor network model.
+
 ** imports
 #+begin_src ipython
   import sys
@@ -887,9 +890,9 @@ from src.analysis.decode import decode_bump
 : 0  2.596293  10.528442 -5.977084        0  0.499  0.0
 : 1  1.005097   4.559610 -5.977084        1  0.499  0.0
 : 2  3.421589  -1.393433 -5.977084        2  0.499  0.0
-: 3  2.504058  -1.709201 -5.977084        3  0.499  0.0
 : 4  2.531308  -2.055951 -5.977084        4  0.499  0.0
 
+: 3  2.504058  -1.709201 -5.977084        3  0.499  0.0
 
 #+begin_src ipython
   res = df.groupby(['time', 'J1'])['rates'].apply(decode_bump).reset_index()