From 532e1ea3438bc795d2f9bf450cd65905146bba6c Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Thu, 3 Oct 2024 11:31:44 +0200
Subject: [PATCH 01/30] Add tutorial
---
doc/tutorial/index.rst | 1 +
doc/tutorial/plasticity.rst | 76 +
python/example/plasticity/03-rates.svg | 2011 +++++++++++++++++++
python/example/plasticity/homeostasis.py | 115 ++
python/example/plasticity/random_network.py | 96 +
python/example/plasticity/unconnected.py | 82 +
6 files changed, 2381 insertions(+)
create mode 100644 doc/tutorial/plasticity.rst
create mode 100644 python/example/plasticity/03-rates.svg
create mode 100644 python/example/plasticity/homeostasis.py
create mode 100644 python/example/plasticity/random_network.py
create mode 100644 python/example/plasticity/unconnected.py
diff --git a/doc/tutorial/index.rst b/doc/tutorial/index.rst
index 8ef283e2d..3adb1567c 100644
--- a/doc/tutorial/index.rst
+++ b/doc/tutorial/index.rst
@@ -78,6 +78,7 @@ Advanced
:maxdepth: 1
nmodl
+ plasticity
Demonstrations
--------------
diff --git a/doc/tutorial/plasticity.rst b/doc/tutorial/plasticity.rst
new file mode 100644
index 000000000..bc9a060d0
--- /dev/null
+++ b/doc/tutorial/plasticity.rst
@@ -0,0 +1,76 @@
+.. _tutorial_plasticity:
+
+In this tutorial, we are going to demonstrate how a network can be built using
+plasticity and homeostatic connection rules. Despite not playing towards Arbor's
+strengths, we choose a LIF (Leaky Integrate and Fire) neuron model, as we are
+primarily interested in examining the required scaffolding.
+
+We will build up the simulation in stages, starting with an unconnected network
+and finishing with a dynamically built connectome.
+
+An Unconnected Network
+----------------------
+
+Consider a collection of ``N`` LIF cells. This will be the starting point for
+our exploration. For now, we set up each cell with a Poissonian input such that
+it will produce spikes periodically at a low frequency.
+
+The Python file ``01-setup.py`` is the scaffolding we will build our simulation
+around and thus contains some passages that might seem redundant now, but will
+be helpful in later steps.
+
+We begin by defining the global settings
+
+.. literalinclude:: ../../python/example/plasticity/unconnected.py
+ :language: python
+ :lines: 9-17
+
+- ``T`` is the total runtime of the simulation in ``ms``
+- ``dT`` defines the _interval_ such that the simulation is advance in discrete
+ steps ``[0, dT, 2 dT, ..., T]``. Later, this will be the timescale of
+ plasticity.
+- ``dt`` is the numerical timestep on which cells evolve
+
+These parameters are used here
+
+.. literalinclude:: ../../python/example/plasticity/step-01.py
+ :language: python
+
+Next, we define the ``recipe`` used to describe the 'network' which is currently
+unconnected.
+
+We also proceed to add spike recording and generating raster/rate plots.
+
+A Randomly Wired Network
+------------------------
+
+We use inheritance to derive a new recipe that contains all the functionality of
+the ``unconnected`` recipe. We add a random connectivity matrix during
+construction, fixed connection weights, and deliver the resulting connections
+via the callback, with the only extra consideration of allowing multiple
+connections between two neurons.
+
+Adding Homeostasis
+------------------
+
+Under the homeostatic model, each cell was a setpoint for the firing rate :math:`\nu^*`
+which is used to determine the creation or destruction of synaptic connections via
+
+.. math::
+
+ \frac{dC}{dt} = \alpha(\nu - \nu^*)
+
+Thus we need to add some extra information to our simulation; namely the
+setpoint :math:`\nu^*_i` for each neuron :math:`i` and the sensitivity parameter
+:math:`\alpha`. We will also use a simplified version of the differential
+equation above, namely adding/deleting exactly one connection if the difference
+of observed to desired spiking frequency exceeds :math:`\pm\alpha`. This is both
+for simplicity and to avoid sudden changes in the network structure.
+
+We do this by tweaking the connection table in between calls to ``run``. In
+particular, we walk the potential pairings of targets and sources in random
+order and check whether the targets requires adding or removing connections. If
+we find an option to fulfill that requirement, we do so and proceed to the next
+target. The randomization is important here, espcially for adding connections as
+to avoid biases, in particular when there are too few eglible connection
+partners.
diff --git a/python/example/plasticity/03-rates.svg b/python/example/plasticity/03-rates.svg
new file mode 100644
index 000000000..c0cbdde28
--- /dev/null
+++ b/python/example/plasticity/03-rates.svg
@@ -0,0 +1,2011 @@
+
+
+
diff --git a/python/example/plasticity/homeostasis.py b/python/example/plasticity/homeostasis.py
new file mode 100644
index 000000000..b296c862d
--- /dev/null
+++ b/python/example/plasticity/homeostasis.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python3
+
+import arbor as A
+from arbor import units as U
+from typing import Any
+import matplotlib.pyplot as plt
+from scipy.signal import savgol_filter
+import numpy as np
+
+from random_network import random_network
+
+# global parameters
+# total runtime [ms]
+T = 10000
+# one interval [ms]
+dT = 100
+# number of intervals
+nT = int((T + dT - 1)//dT)
+# numerical time step [ms]
+dt = 0.1
+
+def randrange(n: int):
+ res = np.arange(n, dtype=int)
+ np.random.shuffle(res)
+ return res
+
+np.random.seed = 23
+
+class homeostatic_network(random_network):
+
+ def __init__(self, N) -> None:
+ super().__init__(N)
+ self.max_inc = 8
+ self.max_out = 8
+ # setpoint rate in kHz
+ self.setpoint = 0.1
+ # sensitivty towards deviations from setpoint
+ self.alpha = 200
+
+
+if __name__ == "__main__":
+ rec = homeostatic_network(10)
+ sim = A.simulation(rec)
+ sim.record(A.spike_recording.all)
+
+ print("Initial network:")
+ print(rec.inc)
+ print(rec.out)
+ print(rec.connections)
+
+ t = 0
+ while t < T:
+ sim.run((t + dT) * U.ms, dt * U.ms)
+ if t < T/2:
+ t += dT
+ continue
+ n = rec.num_cells()
+ rates = np.zeros(n)
+ for (gid, _), time in sim.spikes():
+ if time < t:
+ continue
+ rates[gid] += 1
+ rates /= dT # kHz
+ dC = ((rec.setpoint - rates)*rec.alpha).astype(int)
+ unchangeable = set()
+ added = []
+ deled = []
+ for tgt in randrange(n):
+ if dC[tgt] == 0:
+ continue
+ for src in randrange(n):
+ if dC[tgt] > 0 and rec.add_connection(src, tgt):
+ added.append((src, tgt))
+ break
+ elif dC[tgt] < 0 and rec.del_connection(src, tgt):
+ deled.append((src, tgt))
+ break
+ unchangeable.add(tgt)
+ sim.update(rec)
+ print(f" * t={t:>4} f={rates} [!] {list(unchangeable)} [+] {added} [-] {deled}")
+ t += dT
+
+ print("Final network:")
+ print(rec.inc)
+ print(rec.out)
+ print(rec.connections)
+
+ # Extract spikes
+ times = []
+ gids = []
+ rates = np.zeros(shape=(nT, rec.num_cells()))
+ for (gid, _), time in sim.spikes():
+ times.append(time)
+ gids.append(gid)
+ it = int(time // dT)
+ rates[it, gid] += 1
+
+ fg, ax = plt.subplots()
+ ax.scatter(times, gids, c=gids)
+ ax.set_xlabel('Time $(t/ms)$')
+ ax.set_ylabel('GID')
+ ax.set_xlim(0, T)
+ fg.savefig('03-raster.pdf')
+
+ fg, ax = plt.subplots()
+ ax.plot(np.arange(nT), rates/dT)
+ ax.plot(np.arange(nT), savgol_filter(rates.mean(axis=1)/dT, window_length=5, polyorder=2), color='0.8', lw=4, label='Mean rate')
+ ax.axhline(0.1, label='Setpoint', lw=2, c='0.4')
+ ax.legend()
+ ax.set_xlabel('Interval')
+ ax.set_ylabel('Rate $(kHz)$')
+ ax.set_xlim(0, nT)
+ fg.savefig('03-rates.pdf')
+ fg.savefig('03-rates.png')
+ fg.savefig('03-rates.svg')
diff --git a/python/example/plasticity/random_network.py b/python/example/plasticity/random_network.py
new file mode 100644
index 000000000..921a1598a
--- /dev/null
+++ b/python/example/plasticity/random_network.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+
+import arbor as A
+from arbor import units as U
+from typing import Any
+import matplotlib.pyplot as plt
+import numpy as np
+
+from unconnected import unconnected
+
+# global parameters
+# total runtime [ms]
+T = 100
+# one interval [ms]
+dT = 10
+# number of intervals
+nT = int((T + dT - 1)//dT)
+# numerical time step [ms]
+dt = 0.1
+
+class random_network(unconnected):
+
+ def __init__(self, N) -> None:
+ super().__init__(N)
+ self.syn_weight = 80
+ self.syn_delay = 0.5 * U.ms
+ self.max_inc = 4
+ self.max_out = 4
+ # format [to, from]
+ self.connections = np.zeros(shape=(N, N), dtype=np.uint8)
+ self.inc = np.zeros(N, np.uint8)
+ self.out = np.zeros(N, np.uint8)
+
+ def connections_on(self, gid: int):
+ return [A.connection((source, "src"), "tgt", self.syn_weight, self.syn_delay)
+ for source in range(self.N)
+ for _ in range(self.connections[gid, source])
+ ]
+
+ def add_connection(self, src: int, tgt: int) -> bool:
+ if tgt == src or self.inc[tgt] >= self.max_inc or self.out[src] >= self.max_out:
+ return False
+ self.inc[tgt] += 1
+ self.out[src] += 1
+ self.connections[tgt, src] += 1
+ return True
+
+ def del_connection(self, src: int, tgt: int) -> bool:
+ if tgt == src or self.connections[tgt, src] <= 0:
+ return False
+ self.inc[tgt] -= 1
+ self.out[src] -= 1
+ self.connections[tgt, src] -= 1
+ return True
+
+ def rewire(self):
+ tries = self.N*self.N*self.max_inc*self.max_out
+ while tries > 0 and self.inc.sum() < self.N*self.max_inc and self.out.sum() < self.N*self.max_out:
+ src, tgt = np.random.randint(self.N, size=2, dtype=int)
+ self.add_connection(src, tgt)
+ tries -= 1
+
+
+if __name__ == "__main__":
+ rec = random(10)
+ rec.rewire()
+ sim = A.simulation(rec)
+ sim.record(A.spike_recording.all)
+ t = 0
+ while t < T:
+ t += dT
+ sim.run(t * U.ms, dt * U.ms)
+
+ # Extract spikes
+ times = []
+ gids = []
+ rates = np.zeros(shape=(nT, rec.num_cells()))
+ for (gid, _), time in sim.spikes():
+ times.append(time)
+ gids.append(gid)
+ it = int(time // dT)
+ rates[it, gid] += 1
+
+ fg, ax = plt.subplots()
+ ax.scatter(times, gids, c=gids)
+ ax.set_xlabel('Time $(t/ms)$')
+ ax.set_ylabel('GID')
+ ax.set_xlim(0, T)
+ fg.savefig('02-raster.pdf')
+
+ fg, ax = plt.subplots()
+ ax.plot(np.arange(nT), rates)
+ ax.set_xlabel('Interval')
+ ax.set_ylabel('Rate $(kHz)$')
+ ax.set_xlim(0, nT)
+ fg.savefig('02-rates.pdf')
diff --git a/python/example/plasticity/unconnected.py b/python/example/plasticity/unconnected.py
new file mode 100644
index 000000000..78b848ca7
--- /dev/null
+++ b/python/example/plasticity/unconnected.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+
+import arbor as A
+from arbor import units as U
+from typing import Any
+import matplotlib.pyplot as plt
+import numpy as np
+
+# global parameters
+# total runtime [ms]
+T = 100
+# one interval [ms]
+dT = 10
+# number of intervals
+nT = int((T + dT - 1)//dT)
+# numerical time step [ms]
+dt = 0.1
+
+class unconnected(A.recipe):
+
+ def __init__(self, N) -> None:
+ super().__init__()
+ self.N = N
+ # Cell prototype
+ self.cell = A.lif_cell("src", "tgt")
+ # random seed [0, 100]
+ self.seed = 42
+ # event generator parameters
+ self.gen_weight = 20
+ self.gen_freq = 1 * U.kHz
+
+ def num_cells(self) -> int:
+ return self.N
+
+ def event_generators(self, gid: int) -> list[Any]:
+ return [A.event_generator("tgt",
+ self.gen_weight,
+ A.poisson_schedule(freq=self.gen_freq,
+ seed=self.cell_seed(gid)))]
+
+ def cell_description(self, gid: int) -> Any:
+ return self.cell
+
+ def cell_kind(self, gid: int) -> A.cell_kind:
+ return A.cell_kind.lif
+
+ def cell_seed(self, gid):
+ return self.seed + gid*100
+
+if __name__ == "__main__":
+ rec = unconnected(10)
+ sim = A.simulation(rec)
+ sim.record(A.spike_recording.all)
+
+ t = 0
+ while t < T:
+ t += dT
+ sim.run(t * U.ms, dt * U.ms)
+
+ # Extract spikes
+ times = []
+ gids = []
+ rates = np.zeros(shape=(nT, rec.num_cells()))
+ for (gid, _), time in sim.spikes():
+ times.append(time)
+ gids.append(gid)
+ it = int(time // dT)
+ rates[it, gid] += 1
+
+ fg, ax = plt.subplots()
+ ax.scatter(times, gids, c=gids)
+ ax.set_xlabel('Time $(t/ms)$')
+ ax.set_ylabel('GID')
+ ax.set_xlim(0, T)
+ fg.savefig('01-raster.pdf')
+
+ fg, ax = plt.subplots()
+ ax.plot(np.arange(nT), rates)
+ ax.set_xlabel('Interval')
+ ax.set_ylabel('Rate $(kHz)$')
+ ax.set_xlim(0, nT)
+ fg.savefig('01-rates.pdf')
From dbe7b518e46ccd05cb66b5d86f26c9b444aec115 Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Tue, 8 Oct 2024 10:43:17 +0200
Subject: [PATCH 02/30] Begin code snippets.
---
doc/tutorial/plasticity.rst | 98 +-
python/example/plasticity/01-raster.svg | 1863 ++++++++++++++
python/example/plasticity/01-rates.svg | 2107 ++++++++++++++++
python/example/plasticity/03-rates.svg | 2404 ++++++++++---------
python/example/plasticity/homeostasis.py | 54 +-
python/example/plasticity/random_network.py | 62 +-
python/example/plasticity/unconnected.py | 62 +-
python/example/plasticity/util.py | 45 +
8 files changed, 5404 insertions(+), 1291 deletions(-)
create mode 100644 python/example/plasticity/01-raster.svg
create mode 100644 python/example/plasticity/01-rates.svg
create mode 100644 python/example/plasticity/util.py
diff --git a/doc/tutorial/plasticity.rst b/doc/tutorial/plasticity.rst
index bc9a060d0..d10c079c6 100644
--- a/doc/tutorial/plasticity.rst
+++ b/doc/tutorial/plasticity.rst
@@ -23,32 +23,110 @@ We begin by defining the global settings
.. literalinclude:: ../../python/example/plasticity/unconnected.py
:language: python
- :lines: 9-17
+ :lines: 7-15
+- ``N`` is the cell count of the simulation
- ``T`` is the total runtime of the simulation in ``ms``
-- ``dT`` defines the _interval_ such that the simulation is advance in discrete
- steps ``[0, dT, 2 dT, ..., T]``. Later, this will be the timescale of
+- ``t_interval`` defines the _interval_ such that the simulation is advance in
+ discrete steps ``[0, 1, 2, ...] t_interval``. Later, this will be the timescale of
plasticity.
- ``dt`` is the numerical timestep on which cells evolve
These parameters are used here
-.. literalinclude:: ../../python/example/plasticity/step-01.py
+.. literalinclude:: ../../python/example/plasticity/unconnected.py
+ :language: python
+ :lines: 52-62
+
+where we run the simulation in increments of ``t_interval``.
+
+Back to the recipe; we set a prototypical cell
+
+.. literalinclude:: ../../python/example/plasticity/unconnected.py
:language: python
+ :lines: 23
-Next, we define the ``recipe`` used to describe the 'network' which is currently
-unconnected.
+and deliver it for all ``gid`` s
-We also proceed to add spike recording and generating raster/rate plots.
+.. literalinclude:: ../../python/example/plasticity/unconnected.py
+ :language: python
+ :lines: 42-43
+
+Also, each cell has an event generator attached
+
+.. literalinclude:: ../../python/example/plasticity/unconnected.py
+ :language: python
+ :lines: 33-40
+
+using a Poisson point process seeded with the cell's ``gid``. All other
+parameters are set in the constructor
+
+.. literalinclude:: ../../python/example/plasticity/unconnected.py
+ :language: python
+ :lines: 19-28
+
+We also proceed to add spike recording and generate plots using a helper
+function ``plot_spikes`` from ``util.py``. You can skip the following details
+for now and come back later if you are interested how it works. We generate
+raster plots via ``scatter``. Rates are computed by binning spikes into
+``t_interval`` and the neuron id; the mean rate is the average across the
+neurons smoothed using a Savitzky-Golay filter (``scipy.signal.savgol_filter``).
+We plot per-neuron and mean rates.
A Randomly Wired Network
------------------------
We use inheritance to derive a new recipe that contains all the functionality of
-the ``unconnected`` recipe. We add a random connectivity matrix during
+the ```unconnected`` recipe. We then add a random connectivity matrix during
construction, fixed connection weights, and deliver the resulting connections
-via the callback, with the only extra consideration of allowing multiple
-connections between two neurons.
+via the ``connections_on`` callback, with the only extra consideration of
+allowing multiple connections between two neurons.
+
+In detail, the recipe stores the connection matrix, the current
+incoming/outgoing connections per neuron, and the maximum for both directions
+
+.. literalinclude:: ../../python/example/plasticity/random_network.py
+ :language: python
+ :lines: 26-31
+
+The connection matrix is used to construct connections
+
+.. literalinclude:: ../../python/example/plasticity/random_network.py
+ :language: python
+ :lines: 33-38
+
+together with the fixed connection parameters
+
+.. literalinclude:: ../../python/example/plasticity/random_network.py
+ :language: python
+ :lines: 24-25
+
+We define helper functions ``add|del_connections`` to manipulate the connection
+table while upholding these invariants:
+
+- no self-connections, i.e. ``connection[i, i] == 0``
+- ``inc[i]`` the sum of ``connections[:, i]``
+- no more incoming connections than allowed by ``max_inc``, i.e. ``inc[i] <= max_inc``
+- ``out[i]`` the sum of ``connections[i, :]``
+- no more outgoing connections than allowed by ``max_out``, i.e. ``out[i] <= max_out``
+
+These methods return ``True`` on success and ``False`` otherwise
+
+.. literalinclude:: ../../python/example/plasticity/random_network.py
+ :language: python
+ :lines: 40-54
+
+Both are used in ``rewire`` to produce a random connection matrix
+
+.. literalinclude:: ../../python/example/plasticity/random_network.py
+ :language: python
+ :lines: 56-65
+
+We then proceed to run the simulation and plot the results as before
+
+.. literalinclude:: ../../python/example/plasticity/random_network.py
+ :language: python
+ :lines: 68-79
Adding Homeostasis
------------------
diff --git a/python/example/plasticity/01-raster.svg b/python/example/plasticity/01-raster.svg
new file mode 100644
index 000000000..7b4af213b
--- /dev/null
+++ b/python/example/plasticity/01-raster.svg
@@ -0,0 +1,1863 @@
+
+
+
diff --git a/python/example/plasticity/01-rates.svg b/python/example/plasticity/01-rates.svg
new file mode 100644
index 000000000..67fd888a3
--- /dev/null
+++ b/python/example/plasticity/01-rates.svg
@@ -0,0 +1,2107 @@
+
+
+
diff --git a/python/example/plasticity/03-rates.svg b/python/example/plasticity/03-rates.svg
index c0cbdde28..c6c97fc74 100644
--- a/python/example/plasticity/03-rates.svg
+++ b/python/example/plasticity/03-rates.svg
@@ -6,7 +6,7 @@
- 2024-10-02T20:06:06.006052
+ 2024-10-03T11:36:12.842254
image/svg+xml
@@ -41,12 +41,12 @@ z
-
-
+
@@ -82,7 +82,7 @@ z
-
+
@@ -122,7 +122,7 @@ z
-
+
@@ -157,7 +157,7 @@ z
-
+
@@ -203,7 +203,7 @@ z
-
+
@@ -258,7 +258,7 @@ z
-
+
@@ -445,12 +445,12 @@ z
-
-
+
@@ -475,12 +475,12 @@ z
-
+
-
+
-
+
-
+
@@ -534,12 +534,12 @@ z
-
+
-
+
@@ -550,12 +550,12 @@ z
-
+
-
+
@@ -563,7 +563,89 @@ z
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -678,1143 +760,1143 @@ z
-
-
-
-
-
-
-
-
-
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
+
+
+
+
+
+
+
+
+
-
-
-
+
-
+
-
+
-
-
+
-
+
-
+
+
diff --git a/python/example/plasticity/homeostasis.py b/python/example/plasticity/homeostasis.py
index b296c862d..c15970340 100644
--- a/python/example/plasticity/homeostasis.py
+++ b/python/example/plasticity/homeostasis.py
@@ -4,30 +4,30 @@
from arbor import units as U
from typing import Any
import matplotlib.pyplot as plt
-from scipy.signal import savgol_filter
import numpy as np
+from util import plot_spikes
from random_network import random_network
# global parameters
# total runtime [ms]
T = 10000
# one interval [ms]
-dT = 100
-# number of intervals
-nT = int((T + dT - 1)//dT)
+t_interval = 100
# numerical time step [ms]
dt = 0.1
+
def randrange(n: int):
res = np.arange(n, dtype=int)
np.random.shuffle(res)
return res
+
np.random.seed = 23
-class homeostatic_network(random_network):
+class homeostatic_network(random_network):
def __init__(self, N) -> None:
super().__init__(N)
self.max_inc = 8
@@ -50,9 +50,9 @@ def __init__(self, N) -> None:
t = 0
while t < T:
- sim.run((t + dT) * U.ms, dt * U.ms)
- if t < T/2:
- t += dT
+ sim.run((t + t_interval) * U.ms, dt * U.ms)
+ if t < T / 2:
+ t += t_interval
continue
n = rec.num_cells()
rates = np.zeros(n)
@@ -60,8 +60,8 @@ def __init__(self, N) -> None:
if time < t:
continue
rates[gid] += 1
- rates /= dT # kHz
- dC = ((rec.setpoint - rates)*rec.alpha).astype(int)
+ rates /= t_interval # kHz
+ dC = ((rec.setpoint - rates) * rec.alpha).astype(int)
unchangeable = set()
added = []
deled = []
@@ -78,38 +78,14 @@ def __init__(self, N) -> None:
unchangeable.add(tgt)
sim.update(rec)
print(f" * t={t:>4} f={rates} [!] {list(unchangeable)} [+] {added} [-] {deled}")
- t += dT
+ t += t_interval
print("Final network:")
print(rec.inc)
print(rec.out)
print(rec.connections)
- # Extract spikes
- times = []
- gids = []
- rates = np.zeros(shape=(nT, rec.num_cells()))
- for (gid, _), time in sim.spikes():
- times.append(time)
- gids.append(gid)
- it = int(time // dT)
- rates[it, gid] += 1
-
- fg, ax = plt.subplots()
- ax.scatter(times, gids, c=gids)
- ax.set_xlabel('Time $(t/ms)$')
- ax.set_ylabel('GID')
- ax.set_xlim(0, T)
- fg.savefig('03-raster.pdf')
-
- fg, ax = plt.subplots()
- ax.plot(np.arange(nT), rates/dT)
- ax.plot(np.arange(nT), savgol_filter(rates.mean(axis=1)/dT, window_length=5, polyorder=2), color='0.8', lw=4, label='Mean rate')
- ax.axhline(0.1, label='Setpoint', lw=2, c='0.4')
- ax.legend()
- ax.set_xlabel('Interval')
- ax.set_ylabel('Rate $(kHz)$')
- ax.set_xlim(0, nT)
- fg.savefig('03-rates.pdf')
- fg.savefig('03-rates.png')
- fg.savefig('03-rates.svg')
+ plot_spikes(
+ sim,
+ rec.num_cells(),
+ )
diff --git a/python/example/plasticity/random_network.py b/python/example/plasticity/random_network.py
index 921a1598a..84625afee 100644
--- a/python/example/plasticity/random_network.py
+++ b/python/example/plasticity/random_network.py
@@ -2,40 +2,40 @@
import arbor as A
from arbor import units as U
-from typing import Any
-import matplotlib.pyplot as plt
import numpy as np
+from util import plot_spikes, plot_network
from unconnected import unconnected
# global parameters
+# cell count
+N = 10
# total runtime [ms]
-T = 100
+T = 1000
# one interval [ms]
-dT = 10
-# number of intervals
-nT = int((T + dT - 1)//dT)
+t_interval = 10
# numerical time step [ms]
dt = 0.1
-class random_network(unconnected):
+class random_network(unconnected):
def __init__(self, N) -> None:
super().__init__(N)
self.syn_weight = 80
self.syn_delay = 0.5 * U.ms
- self.max_inc = 4
- self.max_out = 4
# format [to, from]
self.connections = np.zeros(shape=(N, N), dtype=np.uint8)
self.inc = np.zeros(N, np.uint8)
self.out = np.zeros(N, np.uint8)
+ self.max_inc = 4
+ self.max_out = 4
def connections_on(self, gid: int):
- return [A.connection((source, "src"), "tgt", self.syn_weight, self.syn_delay)
- for source in range(self.N)
- for _ in range(self.connections[gid, source])
- ]
+ return [
+ A.connection((source, "src"), "tgt", self.syn_weight, self.syn_delay)
+ for source in range(self.N)
+ for _ in range(self.connections[gid, source])
+ ]
def add_connection(self, src: int, tgt: int) -> bool:
if tgt == src or self.inc[tgt] >= self.max_inc or self.out[src] >= self.max_out:
@@ -54,43 +54,25 @@ def del_connection(self, src: int, tgt: int) -> bool:
return True
def rewire(self):
- tries = self.N*self.N*self.max_inc*self.max_out
- while tries > 0 and self.inc.sum() < self.N*self.max_inc and self.out.sum() < self.N*self.max_out:
+ tries = self.N * self.N * self.max_inc * self.max_out
+ while (
+ tries > 0
+ and self.inc.sum() < self.N * self.max_inc
+ and self.out.sum() < self.N * self.max_out
+ ):
src, tgt = np.random.randint(self.N, size=2, dtype=int)
self.add_connection(src, tgt)
tries -= 1
if __name__ == "__main__":
- rec = random(10)
+ rec = random_network(10)
rec.rewire()
sim = A.simulation(rec)
sim.record(A.spike_recording.all)
t = 0
while t < T:
- t += dT
+ t += t_interval
sim.run(t * U.ms, dt * U.ms)
- # Extract spikes
- times = []
- gids = []
- rates = np.zeros(shape=(nT, rec.num_cells()))
- for (gid, _), time in sim.spikes():
- times.append(time)
- gids.append(gid)
- it = int(time // dT)
- rates[it, gid] += 1
-
- fg, ax = plt.subplots()
- ax.scatter(times, gids, c=gids)
- ax.set_xlabel('Time $(t/ms)$')
- ax.set_ylabel('GID')
- ax.set_xlim(0, T)
- fg.savefig('02-raster.pdf')
-
- fg, ax = plt.subplots()
- ax.plot(np.arange(nT), rates)
- ax.set_xlabel('Interval')
- ax.set_ylabel('Rate $(kHz)$')
- ax.set_xlim(0, nT)
- fg.savefig('02-rates.pdf')
+ plot_spikes(sim, rec.num_cells(), t_interval, T, prefix="02-")
diff --git a/python/example/plasticity/unconnected.py b/python/example/plasticity/unconnected.py
index 78b848ca7..e9a07ec99 100644
--- a/python/example/plasticity/unconnected.py
+++ b/python/example/plasticity/unconnected.py
@@ -2,22 +2,20 @@
import arbor as A
from arbor import units as U
-from typing import Any
-import matplotlib.pyplot as plt
-import numpy as np
+from util import plot_spikes
# global parameters
+# cell count
+N = 10
# total runtime [ms]
-T = 100
+T = 1000
# one interval [ms]
-dT = 10
-# number of intervals
-nT = int((T + dT - 1)//dT)
+t_interval = 10
# numerical time step [ms]
dt = 0.1
-class unconnected(A.recipe):
+class unconnected(A.recipe):
def __init__(self, N) -> None:
super().__init__()
self.N = N
@@ -32,51 +30,33 @@ def __init__(self, N) -> None:
def num_cells(self) -> int:
return self.N
- def event_generators(self, gid: int) -> list[Any]:
- return [A.event_generator("tgt",
- self.gen_weight,
- A.poisson_schedule(freq=self.gen_freq,
- seed=self.cell_seed(gid)))]
+ def event_generators(self, gid: int):
+ return [
+ A.event_generator(
+ "tgt",
+ self.gen_weight,
+ A.poisson_schedule(freq=self.gen_freq, seed=self.cell_seed(gid)),
+ )
+ ]
- def cell_description(self, gid: int) -> Any:
+ def cell_description(self, gid: int):
return self.cell
def cell_kind(self, gid: int) -> A.cell_kind:
return A.cell_kind.lif
- def cell_seed(self, gid):
- return self.seed + gid*100
+ def cell_seed(self, gid: int):
+ return self.seed + gid * 100
+
if __name__ == "__main__":
- rec = unconnected(10)
+ rec = unconnected(N)
sim = A.simulation(rec)
sim.record(A.spike_recording.all)
t = 0
while t < T:
- t += dT
+ t += t_interval
sim.run(t * U.ms, dt * U.ms)
- # Extract spikes
- times = []
- gids = []
- rates = np.zeros(shape=(nT, rec.num_cells()))
- for (gid, _), time in sim.spikes():
- times.append(time)
- gids.append(gid)
- it = int(time // dT)
- rates[it, gid] += 1
-
- fg, ax = plt.subplots()
- ax.scatter(times, gids, c=gids)
- ax.set_xlabel('Time $(t/ms)$')
- ax.set_ylabel('GID')
- ax.set_xlim(0, T)
- fg.savefig('01-raster.pdf')
-
- fg, ax = plt.subplots()
- ax.plot(np.arange(nT), rates)
- ax.set_xlabel('Interval')
- ax.set_ylabel('Rate $(kHz)$')
- ax.set_xlim(0, nT)
- fg.savefig('01-rates.pdf')
+ plot_spikes(sim, rec.num_cells(), t_interval, T, prefix="01-")
diff --git a/python/example/plasticity/util.py b/python/example/plasticity/util.py
new file mode 100644
index 000000000..a089fc626
--- /dev/null
+++ b/python/example/plasticity/util.py
@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import savgol_filter
+
+
+def plot_spikes(sim, n_cells, t_interval, T, prefix=""):
+ # number of intervals
+ n_interval = int((T + t_interval - 1) // t_interval)
+ print(n_interval, T, t_interval)
+
+ # Extract spikes
+ times = []
+ gids = []
+ rates = np.zeros(shape=(n_interval, n_cells))
+ for (gid, _), time in sim.spikes():
+ times.append(time)
+ gids.append(gid)
+ it = int(time // t_interval)
+ rates[it, gid] += 1
+
+ fg, ax = plt.subplots()
+ ax.scatter(times, gids, c=gids)
+ ax.set_xlabel("Time $(t/ms)$")
+ ax.set_ylabel("GID")
+ ax.set_xlim(0, T)
+ fg.savefig(f"{prefix}raster.pdf")
+ fg.savefig(f"{prefix}raster.png")
+ fg.savefig(f"{prefix}raster.svg")
+
+ ts = np.arange(n_interval) * t_interval
+ mean_rate = savgol_filter(
+ rates.mean(axis=1) / t_interval, window_length=5, polyorder=2
+ )
+ fg, ax = plt.subplots()
+ ax.plot(ts, rates)
+ ax.plot(ts, mean_rate, color="0.8", lw=4, label="Mean rate")
+ ax.set_xlabel("Time $(t/ms)$")
+ ax.legend()
+ ax.set_ylabel("Rate $(kHz)$")
+ ax.set_xlim(0, T)
+ fg.savefig(f"{prefix}rates.pdf")
+ fg.savefig(f"{prefix}rates.png")
+ fg.savefig(f"{prefix}rates.svg")
From 8292e816ffeba779c2a337a4800a12042ef566d4 Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Tue, 8 Oct 2024 11:35:24 +0200
Subject: [PATCH 03/30] Add graph plotting.
---
doc/tutorial/plasticity.rst | 85 +-
python/example/plasticity/01-raster.svg | 1592 +-
python/example/plasticity/01-rates.svg | 260 +-
python/example/plasticity/02-graph.svg | 808 +
python/example/plasticity/02-matrix.svg | 388 +
python/example/plasticity/02-raster.svg | 4005 ++++
python/example/plasticity/02-rates.svg | 2090 ++
python/example/plasticity/03-final-graph.svg | 1208 +
python/example/plasticity/03-final-matrix.svg | 388 +
.../example/plasticity/03-initial-graph.svg | 408 +
.../example/plasticity/03-initial-matrix.svg | 388 +
python/example/plasticity/03-raster.svg | 19409 ++++++++++++++++
python/example/plasticity/03-rates.svg | 2982 ++-
python/example/plasticity/homeostasis.py | 52 +-
python/example/plasticity/random_network.py | 3 +-
python/example/plasticity/util.py | 29 +-
16 files changed, 31593 insertions(+), 2502 deletions(-)
create mode 100644 python/example/plasticity/02-graph.svg
create mode 100644 python/example/plasticity/02-matrix.svg
create mode 100644 python/example/plasticity/02-raster.svg
create mode 100644 python/example/plasticity/02-rates.svg
create mode 100644 python/example/plasticity/03-final-graph.svg
create mode 100644 python/example/plasticity/03-final-matrix.svg
create mode 100644 python/example/plasticity/03-initial-graph.svg
create mode 100644 python/example/plasticity/03-initial-matrix.svg
create mode 100644 python/example/plasticity/03-raster.svg
diff --git a/doc/tutorial/plasticity.rst b/doc/tutorial/plasticity.rst
index d10c079c6..7aa3028f1 100644
--- a/doc/tutorial/plasticity.rst
+++ b/doc/tutorial/plasticity.rst
@@ -8,8 +8,18 @@ primarily interested in examining the required scaffolding.
We will build up the simulation in stages, starting with an unconnected network
and finishing with a dynamically built connectome.
-An Unconnected Network
-----------------------
+.. admonition:: Concepts and Requirements
+
+ We cover some advanced topics in this tutorial, mainly structural
+ plasticity. Please refer to other tutorials for the basics of network
+ building. The model employed here --- storing an explicit connection matrix
+ --- is not advisable in most scenarios.
+
+ In addition to Arbor and its requirements, ``scipy``, ``matplotlib``, and
+ ``networkx`` need to be installed.
+
+Unconnected Network
+-------------------
Consider a collection of ``N`` LIF cells. This will be the starting point for
our exploration. For now, we set up each cell with a Poissonian input such that
@@ -128,6 +138,9 @@ We then proceed to run the simulation and plot the results as before
:language: python
:lines: 68-79
+Note that we added a plot of the network connectivity using ``plot_network``
+from ``util`` as well. This generates images of the graph and connection matrix.
+
Adding Homeostasis
------------------
@@ -145,10 +158,64 @@ equation above, namely adding/deleting exactly one connection if the difference
of observed to desired spiking frequency exceeds :math:`\pm\alpha`. This is both
for simplicity and to avoid sudden changes in the network structure.
-We do this by tweaking the connection table in between calls to ``run``. In
-particular, we walk the potential pairings of targets and sources in random
-order and check whether the targets requires adding or removing connections. If
-we find an option to fulfill that requirement, we do so and proceed to the next
-target. The randomization is important here, espcially for adding connections as
-to avoid biases, in particular when there are too few eglible connection
-partners.
+As before, we set up global parameters
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 10-24
+
+and prepare our simulation
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 37-39
+
+Note that our new recipe is almost unaltered from the random network
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 27-33
+
+all changes are contained to the way we run the simulation. To add a further
+interesting feature, we skip the rewiring for the first half of the simulation.
+
+Plasticity is implemented by tweaking the connection table inside the recipe
+between calls to ``run`` and calling ``simulation.update`` with the modified
+recipe:
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 70
+
+Changes are based on the difference of current rate we compute from the spikes
+during the last interval
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 49-54
+
+and the setpoint times the sensitivity
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 55
+
+Then, each potential pairing of target and source is checked in random
+order for whether adding or removing a connection is required
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 59-68
+
+If we find an option to fulfill that requirement, we do so and proceed to the
+next target. The randomization is important here, espcially for adding
+connections as to avoid biases, in particular when there are too few eglible
+connection partners. The ``randrange`` function produces a shuffled range ``[0,
+N)``. We leverage the helper functions from the random network recipe to
+manipulate the connection table, see the discussion above.
+
+Finally, we plot networks and spikes as before
+
+.. literalinclude:: ../../python/example/plasticity/homeostasis.py
+ :language: python
+ :lines: 74-75
diff --git a/python/example/plasticity/01-raster.svg b/python/example/plasticity/01-raster.svg
index 7b4af213b..11b1a2fbd 100644
--- a/python/example/plasticity/01-raster.svg
+++ b/python/example/plasticity/01-raster.svg
@@ -6,7 +6,7 @@
- 2024-10-08T09:16:17.177252
+ 2024-10-08T11:34:50.435338
image/svg+xml
@@ -39,7 +39,7 @@ z
-
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+
-
-
+
+