A Versatile Discrete Simulator

For a case with 3 neurons: N1, N2 and N3, where N1 has N2 and N3 as postsynaptic partners, both with weights 1; N2 has N1 and N3 as postsynaptic partners, both with weights 1; and N3 has N1 has postsynaptic partner with weight 1, the adjacency matrix is:

where we write "." instead of "0" following a convention commonly used for sparse matrix representations. The \(\mathcal{V}_{\cdot \rightarrow i}\) are then the columns of the adjacency matrix, while the \(\mathcal{V}_{i \rightarrow \cdot}\) are its rows. We will also use graphical representations for these matrices as illustrated here:

We can and will generalize the construction of the adjacency matrix by using elements that are not necessarily 0 (no connection) or 1 (connection) but by plugging-in the actual \(w_{i\rightarrow j}\) values.

In the table, the right columns shows the realization \(l_5(i)\) of \(L_5(i)\) for each of the three processes.

Here, 1 is used to indicate the last spike.

Using the previous adjacency matrix, the above realizations lead to:

We see that \(v_{-2}(1)\) is not defined ("?") since \(l_{-2}(1)\) is missing.

This is very "simply" done, we associate to each neuron \(i \in I\) a refractory period \(r_i \in \{1,2,...,r_{\text{max}}\}\) and we replace \eqref{def:potential} by:

then \eqref{def:potential} is a special case of \eqref{def:potential-with-rp} where \(r_i = 1\), \(\forall i \in I\). If we take \(r_i = 2\) for \(i \in \{1,2,3\}\) in Ex. 1.3.1, we get (modified entries appear in red):

The inclusion of a "non trivial" refractctory period (\(r_i > 1\)) in Eq. \eqref{def:potential-with-rp} breaks the Markov property enjoyed by the process \((\mathbb{V}_t)_{t \geq 0 }\) whose dynamics is defined by Eq. \eqref{eq:simple-discrete-V-evolution}. An equivalent of the latter should indeed be written:

We see here that \(V_{t+1}(i)\) depends expicitely on \(L_t(i)\). This means that we must make our state space larger to make room for \(L_t(i)\), namely \((\mathbb{R}^n \times \mathbb{Z}^n)\), in order to get back the Markov property (we need potentially negative values of \(L_t(i)\) to initialize our process, that is why we have \(\mathbb{Z}^n\) and not \(\{0,1,\ldots\}^n\) here). We then have to consider the joint processes: \(\left(\mathbb{V}_t,\mathbb{L}_t\right)_{t \geq 0 }\), where the transitions of \(\mathbb{V}_t\) are given by Eq. \eqref{eq:discrete-V-evolution-with-rp} and the transition of \(\mathbb{L}_t\) are given by:

We can see \(\mathbb{L}_t\) as a vector valued random variable:

This is also "simply" done by adding to each pair \((j,i)\) for which \(w_{j \to i} \neq 0\) a delay \(\tau_{j \to i} \in \{0,1,2,\ldots, \tau_{max}\}\), with \(\tau_{max} < \infty\). Eq. \eqref{def:potential-with-rp} becomes:

If \(r_i = 1\), \(\forall i \in I\), Eq. \eqref{eq:simple-discrete-V-evolution} becomes:

We see that \(V_{t+1 } (i )\) does not depend only on \(\left(V_{t } (1 ), V_{t } (2),\ldots, V_{t } (n )\right)\) anymore, but that we have to go back in time to \(\tau_{max}\). So the nice Markov property we had seems to be lost again.

If we are in the general case, we can consider the joint processes: \[(\mathbb{S}_t)_{t \ge 0} \equiv \left(\mathbb{V}_t,\mathbb{L}_t, \mathbb{X}_t\right)_{t \geq 0 }\, ,\] on the state space \(\mathbb{R}^n \times \mathbb{Z}^n \times \{0,1\}^{n \times (\tau_{max}+1)}\), where \(\mathbb{X}_t\) is defined by:

\(\mathbb{X}_{t}\) is just the subset of all the past events (spike or no spike) generated by the network that can still influence what will happen at \(t+1\).

Then Eq. \eqref{eq:discrete-V-evolution-with-rp} and \eqref{eq:discrete-L-evolution} can be generalized as follows:

We then recover the Markov property.

Notice that the introduction of potentially non-null synaptic delays makes the introduction of self-interaction (\(w_{i \to i} \neq 0\)) meaningful if \(\tau_{i \to i} > r_i\). This can be used to favor a weakly periodic discharge. Imagine for instance that:

• \(\phi_i(v) = \max(0,\min(1,v))\); in other words, \(\phi_i(\,)\) is linear with slope 1 between 0 and 1, is null below 0 and 1 above 1;
• \(r_i = 5\);
• \(w_{i \to i} = 0.7\).

Let us assume that \(i\) spiked at \(t=0\) and that \(i\) is isolated (\(\mathcal{V}_{\cdot{} \to i} = \emptyset\)). Then:

• \(\mathbb{P}\{X_t(i) = 1 \mid L_t(i) \le 5\} = 0\) (refractory period effect);
• \(\mathbb{P}\{X_t(i) = 1 \mid L_t(i) > 5\} = 0.7\).

The probability for observing the next spike at:

• \(t=6\) is 0.7;
• \(t=7\) is \(0.3 \times 0.7\);
• \(t=8\) is \(0.3^2 \times 0.7\);
• …

The distribution of the inter spike interval is therefore a shifted geometric distribution (the shift is 5 and the parameter of the distribution is 0.7). This implies that we have a 0.99 probability for observing the next spike between \(t=6\) and \(t=10\).

If, for all \(i,j \in I\), \(g_{j \to i} ( s) = \rho_i^s\) for all \(s \geq 0\), for some \(\rho_i \in [0,1]\); if \(r_i=1\) and \(\tau_{j \to i} = 0\); then \((\mathbb{V}_t)_{t \geq 0 }\) is a Markov chain whose transitions are given by:

When we will simulate models with leakage that do not correspond to the geometric case of the previous section, we will only consider leakage functions with finite (compact) support. That is, there will be a \(t_{max} < \infty\) such that \(g_{j \to i} (t) = 0\), \(\forall t \ge t_{max}\). Our leakage functions can then be fully specified by vectors in \(\mathbb{R}^{t_{max}}\):

The Markov property is then recovered on the state space: \(\mathbb{R}^n \times \mathbb{Z}^n \times \{0,1\}^{n \times (\tau_{max}+\textcolor{red}{t_{max}}+1)}\) with the process: \((\mathbb{S}_t)_{t \ge 0} \equiv \left(\mathbb{V}_t,\mathbb{L}_t, \mathbb{X}_t\right)_{t \geq 0 }\) where \(\mathbb{X}_t\) is now a bigger matrix than in Eq. \eqref{eq:bigX} since it has to go further back in the past:

A linked list or something close to it, is a uni-dimensional container to which we can add new items to the right (when \(\mathbb{1}_{U_{t+1}(i) \le \phi_i\left(V_t(i)\right)} = 1\)) and from which we can remove an item from the left (\(t-\tau_{max}\) was in the list and we advanced by one time step). This is directly provided by a Python list, the first operation is implemented by the append method and the second by the pop method, as illustrated by the following example:

from random import seed, random
seed(20061001)
tau_max = 5
rho = 0.25
lst = []
t = 0
while t < 20:
    if random() <= rho:
        lst.append(t)
        print("Event at time: " + str(t))
    if len(lst) > 0 and lst[0] < t-tau_max:
        lst.pop(0)
    print("List at time " + str(t) +": " + str(lst))
    t += 1    
    

List at time 0: []
List at time 1: []
List at time 2: []
List at time 3: []
Event at time: 4
List at time 4: [4]
List at time 5: [4]
List at time 6: [4]
List at time 7: [4]
List at time 8: [4]
List at time 9: [4]
List at time 10: []
Event at time: 11
List at time 11: [11]
Event at time: 12
List at time 12: [11, 12]
List at time 13: [11, 12]
List at time 14: [11, 12]
List at time 15: [11, 12]
Event at time: 16
List at time 16: [11, 12, 16]
Event at time: 17
List at time 17: [12, 16, 17]
List at time 18: [16, 17]
List at time 19: [16, 17]

Clearly, a linked list of fixed length can also be used:

from random import seed, random
seed(20061001)
tau_max = 5
rho = 0.25
lst = [0]*(tau_max+1)
t = 0
while t < 20:
    lst.pop(0)
    if random() <= rho:
        lst.append(1)
        print("Event at time: " + str(t))
    else:
        lst.append(0)
    print("List at time " + str(t) +": " + str(lst))
    t += 1    
    

List at time 0: [0, 0, 0, 0, 0, 0]
List at time 1: [0, 0, 0, 0, 0, 0]
List at time 2: [0, 0, 0, 0, 0, 0]
List at time 3: [0, 0, 0, 0, 0, 0]
Event at time: 4
List at time 4: [0, 0, 0, 0, 0, 1]
List at time 5: [0, 0, 0, 0, 1, 0]
List at time 6: [0, 0, 0, 1, 0, 0]
List at time 7: [0, 0, 1, 0, 0, 0]
List at time 8: [0, 1, 0, 0, 0, 0]
List at time 9: [1, 0, 0, 0, 0, 0]
List at time 10: [0, 0, 0, 0, 0, 0]
Event at time: 11
List at time 11: [0, 0, 0, 0, 0, 1]
Event at time: 12
List at time 12: [0, 0, 0, 0, 1, 1]
List at time 13: [0, 0, 0, 1, 1, 0]
List at time 14: [0, 0, 1, 1, 0, 0]
List at time 15: [0, 1, 1, 0, 0, 0]
Event at time: 16
List at time 16: [1, 1, 0, 0, 0, 1]
Event at time: 17
List at time 17: [1, 0, 0, 0, 1, 1]
List at time 18: [0, 0, 0, 1, 1, 0]
List at time 19: [0, 0, 1, 1, 0, 0]

The array object of the array module also provide the necessary methods append and pop needed to reproduce the last exemple, so lets try:

from array import array
from random import seed, random
seed(20061001)
tau_max = 5
rho = 0.25
lst = array('B',[0]*(tau_max+1))
t = 0
while t < 20:
    lst.pop(0)
    if random() <= rho:
        lst.append(1)
        print("Event at time: " + str(t))
    else:
        lst.append(0)
    print("List at time " + str(t) +": " + str(list(lst)))
    t += 1    
    

List at time 0: [0, 0, 0, 0, 0, 0]
List at time 1: [0, 0, 0, 0, 0, 0]
List at time 2: [0, 0, 0, 0, 0, 0]
List at time 3: [0, 0, 0, 0, 0, 0]
Event at time: 4
List at time 4: [0, 0, 0, 0, 0, 1]
List at time 5: [0, 0, 0, 0, 1, 0]
List at time 6: [0, 0, 0, 1, 0, 0]
List at time 7: [0, 0, 1, 0, 0, 0]
List at time 8: [0, 1, 0, 0, 0, 0]
List at time 9: [1, 0, 0, 0, 0, 0]
List at time 10: [0, 0, 0, 0, 0, 0]
Event at time: 11
List at time 11: [0, 0, 0, 0, 0, 1]
Event at time: 12
List at time 12: [0, 0, 0, 0, 1, 1]
List at time 13: [0, 0, 0, 1, 1, 0]
List at time 14: [0, 0, 1, 1, 0, 0]
List at time 15: [0, 1, 1, 0, 0, 0]
Event at time: 16
List at time 16: [1, 1, 0, 0, 0, 1]
Event at time: 17
List at time 17: [1, 0, 0, 0, 1, 1]
List at time 18: [0, 0, 0, 1, 1, 0]
List at time 19: [0, 0, 1, 1, 0, 0]

We can choose to work with a fixed length list on which we implement a circular buffer as shown on stackoverflow that we modify as follows.

class circularlist(object):
    from array import array
    def __init__(self, size):
        """Initialization"""
        self.index = 0
        self.size = size
        self._data = array('B',[0]*size)

    def append(self, value):
        """Append an element"""
        self._data[self.index] = value
        self.index = (self.index + 1) % self.size

    def __getitem__(self, key):
        """Get element by index, relative to the current index"""
        return(self._data[(key + self.index) % self.size])
        
    def __repr__(self):
        """Return string representation"""
        return str([self[i] for i in range(self.size)])

With this class defined, the previous example becomes:

from random import seed, random
seed(20061001)
tau_max = 5
rho = 0.25
buffer = circularlist(tau_max+1)
t = 0
while t < 20:
    if random() <= rho:
        print("Event at time: " + str(t))
        buffer.append(1)
    else:
        buffer.append(0)
    print("Buffer at time " + str(t) +": " + str(buffer))
    t += 1    
    

Buffer at time 0: [0, 0, 0, 0, 0, 0]
Buffer at time 1: [0, 0, 0, 0, 0, 0]
Buffer at time 2: [0, 0, 0, 0, 0, 0]
Buffer at time 3: [0, 0, 0, 0, 0, 0]
Event at time: 4
Buffer at time 4: [0, 0, 0, 0, 0, 1]
Buffer at time 5: [0, 0, 0, 0, 1, 0]
Buffer at time 6: [0, 0, 0, 1, 0, 0]
Buffer at time 7: [0, 0, 1, 0, 0, 0]
Buffer at time 8: [0, 1, 0, 0, 0, 0]
Buffer at time 9: [1, 0, 0, 0, 0, 0]
Buffer at time 10: [0, 0, 0, 0, 0, 0]
Event at time: 11
Buffer at time 11: [0, 0, 0, 0, 0, 1]
Event at time: 12
Buffer at time 12: [0, 0, 0, 0, 1, 1]
Buffer at time 13: [0, 0, 0, 1, 1, 0]
Buffer at time 14: [0, 0, 1, 1, 0, 0]
Buffer at time 15: [0, 1, 1, 0, 0, 0]
Event at time: 16
Buffer at time 16: [1, 1, 0, 0, 0, 1]
Event at time: 17
Buffer at time 17: [1, 0, 0, 0, 1, 1]
Buffer at time 18: [0, 0, 0, 1, 1, 0]
Buffer at time 19: [0, 0, 1, 1, 0, 0]

We consider here a genuine binary coding of the rows of \(\mathbf{x}_t\), a realisation of \(\mathbb{X}_t\). More precisely, the above example at time 17 with a \(\tau_{max} = 5\) corresponds to: \[\left(x_{12}(i)=1, x_{13}(i)=0, x_{14}(i)=0, x_{15}(i)=0, x_{16}(i)=1, x_{17}(i)=1\right)\] We can represent the binary sequence 100011 by the integer 35: \[35 = 1 \times 2^5 + 0 \times 2^4 + 0 \times 2^3 + 0 \times 2^2 + 1 \times 2^1 + 1 \times 2^0\, ,\] in Python the conversion is automatically performed by prefixing the sequence with "0b":

x_17 = 0b100011
x_17

35

From this point, advancing by one time step is simply done by multiplying the number by 2 (a left shift 35 << 1) before taking the modulus with \(2^6\):

x_18 = (x_17 << 1) % 2**6
print("{0:0>6b}".format(x_18))

000110

If a spike occurs, we just have to add 1. Then when we write:

t = 17 
seq = 0b110001
seq <<= 1
seq %= 2**6
seq += 1
print("Sequence at time {0:2d} : {1:0>6b}".format(t,seq))

Sequence at time 17 : 100011

the sequence "100011" should be read from right to left. Time is 17 and the rightmost element is 1, so there was a spike at time 17. The second rightmost element is 1, so there was also a spike at time 17-1=16, the next three rightmost elements are 0, so there were no spike at times 17-2=15, 17-3=14 or 17-4=13; the leftmost element is 1, so there was a spike at time 17-5=12. In order to know quickly if there was a spike s time step before t (s=0,1,2,3,4,5) we can simply use a binary left shift >> operation with t-s steps and check if the remainder of a division by 2 is 1:

t = 17 
print("The present time is " + str(t))
print("Was there a spike at time " + str(t) + "?")
print("  " + str((x_17 >> (t-t)) % 2 == 1))
print("Was there a spike at time " + str(t-1)+ "?")
print("  " + str((x_17 >> (t-t+1)) % 2 == 1))
print("Was there a spike at time " + str(t-2)+ "?")
print("  " + str((x_17 >> (t-t+2)) % 2 == 1))

The present time is 17
Was there a spike at time 17?
  True
Was there a spike at time 16?
  True
Was there a spike at time 15?
  False

So let us re-code our previous example:

from random import seed, random
seed(20061001)
tau_max = 5
nb_seq = 2**(tau_max+1)
rho = 0.25
seq = 0b0
t = 0
while t < 20:
   seq <<= 1
   seq %= nb_seq
   if random() <= rho:
      seq += 1
      print("Event at time       : " + str(t))
   print("Sequence at time {0:2d} : {1:0>6b}".format(t,seq))
   t += 1    
   

Sequence at time  0 : 000000
Sequence at time  1 : 000000
Sequence at time  2 : 000000
Sequence at time  3 : 000000
Event at time       : 4
Sequence at time  4 : 000001
Sequence at time  5 : 000010
Sequence at time  6 : 000100
Sequence at time  7 : 001000
Sequence at time  8 : 010000
Sequence at time  9 : 100000
Sequence at time 10 : 000000
Event at time       : 11
Sequence at time 11 : 000001
Event at time       : 12
Sequence at time 12 : 000011
Sequence at time 13 : 000110
Sequence at time 14 : 001100
Sequence at time 15 : 011000
Event at time       : 16
Sequence at time 16 : 110001
Event at time       : 17
Sequence at time 17 : 100011
Sequence at time 18 : 000110
Sequence at time 19 : 001100

We can try doing the same thing using the numpy types and functions:

import numpy as np
rng = np.random.default_rng(20061001)
tau_max = 5
nb_seq = 2**(tau_max+1)
rho = np.float64(0.25)
seq = np.uint8(0)
n_past = np.uint8(2**(tau_max+1)-1)
one = np.uint8(1)
t = 0
while t < 20:
   seq = np.left_shift(seq,one)
   seq = np.bitwise_and(seq,n_past)
   if rng.random()  <= rho:
      seq += 1
      print("Event at time       : " + str(t))
   print("Sequence at time {0:2d} : {1:0>6b}".format(t,seq))
   t += 1    
   

Sequence at time  0 : 000000
Event at time       : 1
Sequence at time  1 : 000001
Sequence at time  2 : 000010
Sequence at time  3 : 000100
Event at time       : 4
Sequence at time  4 : 001001
Sequence at time  5 : 010010
Event at time       : 6
Sequence at time  6 : 100101
Event at time       : 7
Sequence at time  7 : 001011
Sequence at time  8 : 010110
Sequence at time  9 : 101100
Event at time       : 10
Sequence at time 10 : 011001
Sequence at time 11 : 110010
Sequence at time 12 : 100100
Sequence at time 13 : 001000
Sequence at time 14 : 010000
Sequence at time 15 : 100000
Sequence at time 16 : 000000
Sequence at time 17 : 000000
Sequence at time 18 : 000000
Sequence at time 19 : 000000

Take \(\tau_{max} = 50\), \(2^{20}\) steps, a success probability of 0.01 (making IID draws as usual). Initialize your "past" as an empty list when using a variable length scheme (first case of Sec. 2.4.1) or with zeros for a fixed length scheme. At each time step draw a result, make the "past" evolve and draws to integers, j and k, uniformly from \(\{0,1,\ldots, \tau_{max}\}\), before checking if the past contains a spike at time \(t-j\) and \(t-k\). Consider all the above schemes as alternatives (always starting from the same seed) and use timeit to compare them. As an example for the fixed length linked list described in Sec. 2.4.1, you should write something like:

def sim_with_fixed_length_linked_list(n_steps=2**20,
                                      rho_success=0.01,
                                      tau_max=50):
    from random import random, randint
    past = [0]*(tau_max+1)
    t = 0
    found_in_past = 0
    while t < n_steps:
        past.pop(0)
        if random() <= rho_success:
            past.append(1)
        else:
            past.append(0)
        j = randint(0,tau_max)
        found_in_past += past[tau_max-j]
        k = randint(0,tau_max)
        found_in_past += past[tau_max-k]
        t += 1
    print("Found: " + str(found_in_past))    
    return found_in_past

import timeit
from random import seed
seed(20061001)
print(timeit.timeit('sim_with_fixed_length_linked_list()', number=10, globals=globals()))

Found: 20773
Found: 21312
Found: 21316
Found: 20974
Found: 21268
Found: 21059
Found: 20887
Found: 20812
Found: 21493
Found: 20955
18.371022567007458

The code prints the number of spike founds while exploring the past (this is done twice after each time step, making an expected value for this number of \(2 \times 0.01 \times 2^{20} = 20971.52\)). This is a way to (partially) check that the different implementations do generate exactly the same sequences (as long as they start from the same seed).

We are considering here \(2^{20}\) IID draws from a binomial distribution with a parameter of 0.01, we therefore expect:

2**20*0.01

10485.76

successes with a variance of:

2**20*0.01*(1-0.01)

10380.9024

The z-score is:

from math import sqrt
(sum(ref_lst) - 2**20*0.01)/sqrt(2**20*0.01*(1-0.01))

-0.7141264788892625

Our simulation is compatible with a draw from a Binomial distribution with \(n=2^{20}\) and \(p=0.01\).

We simulate idx_lst as follows:

from random import seed, random
import sys
seed(20061001)
idx_lst = [i for i in range(2**20) if random() <= 0.01]
print(total_size(idx_lst))

376740

We compare the contents with:

all([ref_lst[i] for i in idx_lst])

True

The size comparison with getsizeof:

sys.getsizeof(ref_lst)/sys.getsizeof(idx_lst)

99.19141542218465

This ratio is what we expect considering that only a fraction of 0.01 elements of ref_lst are equal to 1. Therefore, the memory size of idx_lst is very close to the size of an array of pointers (with an 8 bytes addressable space) with len(idx_lst) elements:

sys.getsizeof(idx_lst)/(8*len(idx_lst))

1.0224719101123596

If we compare the real sizes we get:

total_size(ref_lst)/total_size(idx_lst)

22.42602325211021

The ratio is smaller because the objects to which the elements of the idx_lst list are pointing to are all different integers leading to a real memory size of: len(idx_lst) * (8+28) since the memory taken by an integer in Python (at least on my machine) is 28 Bytes. 28 Bytes is really a lot for an integer; it is propbably because there is no size limit for integers in Python.

We store ref_lst in an array of array with an unsigned char encoding:

from array import array
ref_array = array('B',ref_lst)
total_size(ref_array)

1048656

This is indeed what we expect from \(2^{20}\) elements taking 1 Byte each. The ratio with our reference is:

total_size(ref_lst)/total_size(ref_array)

8.056769808211653

Again that makes sense: every element of ref_lst is a pointer (taking 8 Bytes) to one of 2 immutable objects (True or False), while every element of ref_array takes 1 Byte.

We now store idx_lst in an array of array with an unsigned int encoding:

idx_array = array('I',idx_lst)
total_size(idx_array)

41732

The ratio with our reference is:

total_size(ref_lst)/total_size(idx_array)

202.45327326751652

That's getting huge!

We store ref_lst in a Boolean numpy array with:

import numpy as np
ref_np = np.array(ref_lst,dtype=np.bool_)
total_size(ref_np)

1048688

You can check that getsizeof and total_size give here the same results (as epxpected, since numpy arrays are not arrays of pointers).

We see that ref_np takes very nearly 1 Byte per element:

total_size(ref_np)/2**20

1.0001068115234375

The size comparison with ref_lst gives:

total_size(ref_lst)/total_size(ref_np)

8.056523961368873

Same as what we get using an array from the array module and an unsigned char encoding.

We store ref_lst in numpy packbits array with:

ref_packed = np.packbits(ref_lst)
total_size(ref_packed)

131184

This is expected since each element of ref_packed takes 1 Byte of memory and correspond to 8 successive binary/Boolean elements of ref_lst, the memory size should therefore be very close:

2**20//8

131072

The size comparison with ref_lst gives:

total_size(ref_lst)/total_size(ref_packed)

64.40404317599707

We now that 8 successive times are included in a single element of ref_packed, we should therefore get first the result, k, of the target time \(t\) divided (integer division) by 8 (k = t//8), get element k of ref_packed and check if the latter is bit 7-(t-k*8) is one (check stackoverflow Get n-th bit of an integer to see how to do that if you don't know). For the first actual spike time we get:

t1 = idx_lst[0]
k1 = t1 // 8
r1 = t1 - k1*8
(ref_packed[k1] >> 7-r1) == 1

True

For the time immediately before:

t2 = idx_lst[0]-1
k2 = t2 // 8
r2 = t2 - k2*8
(ref_packed[k2]  >> 7-r2) == 1

False

As an alternative, you can use unpackbits as follows:

(np.unpackbits(ref_packed[k1])[r1] == 1, np.unpackbits(ref_packed[k2])[r2] == 1) 

(True, False)

We a function run_sim_iid doing the whole simulation plus writing as follows:

def run_sim_iid(rho, ntw_size, n_steps, seed, fname):
    """Runs an IID simulation of a network with 'ntw_size' neurons
    a success probability 'rho' for 'n_steps' seeding the PRNG to
    'seed' and storing the results in file 'fname'

    Parameters
    ----------
    rho: a real between 0 and 1
    ntw_size: a positive integer
    n_steps: a positive integer
    seed: the seed
    fname: a string

    Returns
    -------
    Nothing, the function is used for its side effect: a file is opened
    and the simulation results are stored in it.
    """
    from random import seed
    if not isinstance(n_steps,int):
        raise TypeError("n_steps should be an integer.")
    if not n_steps > 0:
        raise ValueError("Expecting: n_steps > 0.")
    sgl_step_iid = mk_sgl_step_iid(rho,ntw_size)
    with open(fname,"w") as f:
        for i in range(n_steps):
            f.write(' '.join(map(str, sgl_step_iid()))+'\n')

We can now time it in the same conditions as before:

print(timeit.timeit('run_sim_iid(0.01, 10**3, 2**20, 20110928, "sim_iid_t1")',
                    number=1, globals=globals()))

66.43105684299371

Two extra seconds is not much to pay.

We a function run_sim_iid_np doing the whole simulation plus writing as follows:

def run_sim_iid_np(rho, ntw_size, n_steps, seed, fname):
    """Runs an IID simulation of a network with 'ntw_size' neurons
    a success probability 'rho' for 'n_steps' seeding the PRNG to
    'seed' and storing the results in file 'fname'

    Parameters
    ----------
    rho: a real between 0 and 1
    ntw_size: a positive integer
    n_steps: a positive integer
    seed: the seed
    fname: a string

    Returns
    -------
    Nothing, the function is used for its side effect: a file is opened
    and the simulation results are stored in it.
    """
    if not isinstance(n_steps,int):
        raise TypeError("n_steps should be an integer.")
    if not n_steps > 0:
        raise ValueError("Expecting: n_steps > 0.")
    sgl_step_iid_np = mk_sgl_step_iid_np(rho,ntw_size,seed)
    with open(fname,"w") as f:
        for i in range(n_steps):
            f.write(' '.join(map(str, sgl_step_iid_np()))+'\n')

We can now time it in the same conditions as before:

print(timeit.timeit('run_sim_iid_np(0.01, 10**3, 2**20, 20110928, "sim_iid_t2")',
                    number=1, globals=globals()))

15.068215111008612

Again, two extra seconds is not much to pay. Keep in mind that we have a very fast "data" generation in this latter case, while in the real setting we will simulate we are going to be closer to the first.

The major "subtlety" is finding out if the indices are present or not, this is done with the in operator:

def sim_with_variable_length_linked_list(n_steps=2**20,
                                         rho_success=0.01,
                                         tau_max=50):
    from random import random, randint
    past = []
    t = 0
    found_in_past = 0
    while t < n_steps:
        if random() <= rho_success:
            past.append(t)
        if len(past) > 0 and past[0] < t-tau_max:
            past.pop(0)
        j = randint(0,tau_max)
        found_in_past += int(t - j in past)
        k = randint(0,tau_max)
        found_in_past += int(t - k in past)
        t += 1
    print("Found: " + str(found_in_past))    
    return found_in_past

seed(20061001)
print(timeit.timeit('sim_with_variable_length_linked_list()', number=10, globals=globals()))

Found: 20773
Found: 21312
Found: 21316
Found: 20974
Found: 21268
Found: 21059
Found: 20887
Found: 20812
Found: 21493
Found: 20955
20.48626351598068

def sim_with_array_based_fixed_length_linked_list(n_steps=2**20,
                                                  rho_success=0.01,
                                                  tau_max=50):
    from random import random, randint
    from array import array
    past = array('B',[0]*(tau_max+1))
    t = 0
    found_in_past = 0
    while t < n_steps:
        past.pop(0)
        if random() <= rho_success:
            past.append(1)
        else:
            past.append(0)
        j = randint(0,tau_max)
        found_in_past += past[tau_max-j]
        k = randint(0,tau_max)
        found_in_past += past[tau_max-k]
        t += 1
    print("Found: " + str(found_in_past))    
    return found_in_past

import timeit
from random import seed
seed(20061001)
print(timeit.timeit('sim_with_array_based_fixed_length_linked_list()',
                    number=10, globals=globals()))

Found: 20773
Found: 21312
Found: 21316
Found: 20974
Found: 21268
Found: 21059
Found: 20887
Found: 20812
Found: 21493
Found: 20955
18.992268898990005

def sim_with_circularlist(n_steps=2**20,
                          rho_success=0.01,
                          tau_max=50):
    from random import random, randint
    past = circularlist(tau_max+1)
    t = 0
    found_in_past = 0
    while t < n_steps:
        if random() <= rho_success:
            past.append(1)
        else:
            past.append(0)
        j = randint(0,tau_max)
        found_in_past += past[tau_max-j]
        k = randint(0,tau_max)
        found_in_past += past[tau_max-k]
        t += 1
    print("Found: " + str(found_in_past))    
    return found_in_past

seed(20061001)
print(timeit.timeit('sim_with_circularlist()', number=10, globals=globals()))

Found: 20773
Found: 21312
Found: 21316
Found: 20974
Found: 21268
Found: 21059
Found: 20887
Found: 20812
Found: 21493
Found: 20955
23.96578602300724

def sim_with_binary_past(n_steps=2**20,
                         rho_success=0.01,
                         tau_max=50):
    from random import random, randint
    past = 0
    n_past = 2**(tau_max+1)
    t = 0
    found_in_past = 0
    while t < n_steps:
        past <<= 1
        past %= n_past
        if random() <= rho_success:
            past += 1
        j = randint(0,tau_max)
        found_in_past += (past >> j) % 2
        k = randint(0,tau_max)
        found_in_past += (past >> k) % 2
        t += 1
    print("Found: " + str(found_in_past))    
    return found_in_past

seed(20061001)
print(timeit.timeit('sim_with_binary_past()', number=10, globals=globals()))

Found: 20773
Found: 21312
Found: 21316
Found: 20974
Found: 21268
Found: 21059
Found: 20887
Found: 20812
Found: 21493
Found: 20955
19.031911543017486

We see that the best solutions with this parameter set are the fixed length linked list and the binary coding.

For the records, a numpy version, be careful, the function is called only once by the timer!:

def sim_with_binary_past_np(n_steps=2**20,
                            rho_success=0.01,
                            tau_max=50,
                            seed=20061001):
    import numpy as np
    rho_success = np.float64(rho_success)
    rng = np.random.default_rng(seed)
    past = np.uint64(0)
    n_past = np.uint64(2**(tau_max+1)-1)
    one = np.uint64(1)
    t = 0
    found_in_past = 0
    while t < n_steps:
        past = np.left_shift(past,one)
        past = np.bitwise_and(past,n_past)
        if rng.random() <= rho_success:
            past += one
        j = rng.integers(0,tau_max,dtype=np.uint64,endpoint=True)
        found_in_past += int(np.bitwise_and(np.right_shift(past, j),one))
        k = rng.integers(0,tau_max,dtype=np.uint64,endpoint=True)
        found_in_past += int(np.bitwise_and(np.right_shift(past, k),one))
        t += 1
    print("Found: " + str(found_in_past))    
    return found_in_past

import timeit
print(timeit.timeit('sim_with_binary_past_np()', number=1, globals=globals()))

Found: 21591
14.47421620201203

Going back to Sec. 1.5 and 1.6, we see that in the general case a model specification must include:

n

the number of neurons in the network (a positive integer).

edges/synapses properties

\(\mathbb{W}\)

an n x n matrix of synaptic weights (real elements).

\(\tau_{max}\)

the largest synaptic delay (a positive or null integer).

\((\tau_{i \to j})_{i,j \in \{1,\ldots,n\}}\)

an n x n matrix of synaptic delays (positive of null integer elements).

nodes/neurons properties

\(r_{max}\)

the largest refractory period (a positive integer).

\((r_i)_{i \in \{1,\ldots,n\}}\)

a vector of refractory periods (n positive integer elements).

\((\phi_i)_{i \in \{1,\ldots,n\}}\)

a vector of n spiking probability functions.

\((g_{i \to j})_{i,j \in \{1,\ldots,n\}}\)

an n x n matrix a synaptic leakage functions

In order to advance in time our simulation code must keep track, in the most general case, of (Sec. 1.5.2):

\(\mathbf{x}_{t}\)

the realization of the "window on the past", \(\mathbb{X}_{t}\), defined by Eq. \eqref{eq:bigX}, a \(n \times \tau_{max}\) matrix of a vector/list of linked-lists / circular buffers / binary encodings (Sec. 2.4).

\(\mathbf{v}_{t}\)

the realization of the present membrane potentials, \(\mathbb{V}_{t}\), defined by Eq. \eqref{eq:bigV}, a vector/list with \(n\) elements.

\(\mathbf{l}_{t}\)

the realization of the present "elapsed times since the last spike", \(\mathbb{L}_{t}\), defined by Eq. \eqref{eq:bigL}, a vector/list with \(n\) elements.

Notice that all our "evolution equations"—Eq. \eqref{eq:simple-discrete-V-evolution}, \eqref{def:potential-with-rp}, \eqref{eq:complex-discrete-V-evolution}, \eqref{eq:discrete-V-at-t+1-given}, \eqref{def:potentialwithg}—, consider each neuron \(i\) of the network in turn and, for each \(i\), involve a summation over all the neurons of network. This latter summation can nevertheless be made more efficient if we restrict it to the presynaptic neurons to neuron \(i\): the ones for which \(w_{j \to i} \neq 0\) (or the elements of the set \({\mathcal V}_{ \cdot \to i }\) defined in Sec. 1.1). It is therefore a good idea to include in our model specification:

\(\left({\mathcal V}_{ \cdot \to i }\right)_{i \in \{1,\ldots,n\}}\)

a vector/list of lists.

The simulation code will need to have access to these long lists of elements (of very different types). This means that the functions we will have to define will need a lot of input parameters. One way to make our code tidier is then to define a Class whose attributes are the elements we just listed as well as private functions, called methods in the object oriented programming paradigm. An understanding of the class/method mechanism of Python, at the level of Chap. 2 of Goodrich, Tamassia and Goldwasser, Data Structures and Algorithm in Python (2013)—that is, slightly above the level of the Classes section of the official tutorial—, is required in what follows.

We do not necessarily want to have as many different \(\phi_i\) functions as we have neurons in the network or as many different \(g_{i \to j}\) functions as we have actual synapses. In real networks, neurons come in types and in our highly simplified representation it makes sense to identify a type with a rate function, like a \(\phi_{E}\) for excitatory neurons and a \(\phi_{I}\) for inhibitory ones (the latter might be refined if we wanted to add diversity in the inhibitory neurons population). A type would also be associated with a refractory period (\(r_E\) and \(r_i\)) and a self-interaction (if we decide to work with one, \(g_{E}\) and \(g_{I}\)). We would also have as many types of coupling functions \(g_{j \to i}\) for \(j \ne i\) as we have pairs of types: \(E\to E, E \to I, I \to I, I \to I\).

We can therefore reduce the memory used to store a model instance by attaching a label to each neuron specifying its type or kind from a finite set, like \(\mathcal{K} \equiv \{E,I\}\). We then need as many refractory periods and self-interaction functions as we have types and as many coupling functions as the number of types squared. Think about it if we have a network made of 10000 neurons made of 2 types, the potential (memory) saving is huge. We can therefore make our previous "nodes/neurons properties" list more specific by using the types/kinds:

types

a tuple representing the \(\mathcal{K}\) set above (e.g. ("E","I")).

\((r_k)_{k \in \mathcal{K}}\)

the type specific refractory periods.

\(r_{max}\)

\(r_{max} \equiv \max_{k \in \mathcal{K}} r_k\).

\((\phi_k)_{k \in \mathcal{K}}\)

a vector of \(|\mathcal{K}|\) (the cardinality of set \(\mathcal{K}\)) spiking probability functions.

\((g_k)_{k \in \mathcal{K}}\)

a vector of \(|\mathcal{K}|\) self-interaction functions.

\((g_{i \to j})_{i,j \in \mathcal{K}}\)

a \(|\mathcal{K}| \times |\mathcal{K}|\) matrix of coupling functions.

Let us for instance consider a simple \(\phi(v)\) function that is linear between \(v_{min}\) and \(v_{max}\) going from 0 to 1, 0 below \(v_{min}\) and 1 above \(v_{max}\). A simple implementation in Python would be:

def phi_linear(v, v_min, v_max):
    if v < v_min: return 0
    if v > v_max: return 1
    return (v-v_min)/(v_max-v_min)

Notice that we deliberately avoid parameters checking (v_max should be > v_min for instance) in order to gain speed. This function could be replaced by:

def phi_linear2(v, v_min, v_max):
    return min(1,max(0,(v-v_min)/(v_max-v_min)))

These are two different ways to define in Python the same mathematical function.

Compare the times required to evaluate these two functions, as well as a direct evaluation of min(1,max(0,(v-v_min)/(v_max-v_min))) with v_min = 0, v_max = 1 and v taking values in [i/25 fir i in range(26).

In the following questions we will also make the comparison with a slightly more costly \(\phi(\,)\) involving an hyperbolic tangent:

def phi_tanh(v, offset=5, slope=1):
    from math import tanh
    return (1.0+tanh((v-offset)*slope))/2

It is possible and easy to define C functions that do the same jobs as phi_linear and phi_tanh, to compile them (if you have a C compiler like gcc) and to call them from Python thanks to the functionalities of the ctypes module (this could also be done, resulting in principle in faster codes, with SWIG or Cython, but that's heavy, for me). We could for instance define function phi_c as follows and store this definition in file phi_c.c:

#include <stdlib.h>
#include <math.h>

double phi_l(double v, double v_min, double v_max) {
  if (v < v_min) return 0.0;
  if (v > v_max) return 1.0;
  return (v-v_min)/(v_max-v_min);
}

double phi_t(double v, double offset, double slope) {
  return (1.0+tanh((v-offset)*slope))*0.5;
}    

We compile the function into a dynamic library with:

gcc -O2 -shared -o phi_c.so -fPIC phi_c.c

In our Python session we import ctypes and 'load' the shared library:

import ctypes as ct
_phi_clib = ct.cdll.LoadLibrary('./phi_c.so')

We now have to look at the signature of the functions we want to call, like double phi_c(double v, double v_min, double v_max). Our functions takes 3 floating point numbers in double precision and returns one floating point number is double precision. If the equivalent Python types are not specified, Python assumes essentially that everything is an integer and that's neither the case for the parameters nor for the returned value. We therefore do:

_phi_clib.phi_l.argtypes=[ct.c_double,ct.c_double,ct.c_double]  # Parameter types
_phi_clib.phi_l.restype = ct.c_double            # Returned value type
_phi_clib.phi_t.argtypes=[ct.c_double,ct.c_double,ct.c_double]  # Parameter types
_phi_clib.phi_t.restype = ct.c_double            # Returned value type

We can now check that _philib.phi_l gives the same values as phi_linear:

all([_phi_clib.phi_l(i/25,0,1) == phi_linear(i/25,0,1) for i in range(26)])

True

And that _philib.phi_t gives the same values as phi_tanh:

all([_phi_clib.phi_t(i/5,0,1) == phi_tanh(i/5,0,1) for i in range(26)])

True

Make a time comparison between phi_linear and _phi_clib.phi_l, then between phi_tanh and _phi_clib.phi_t like in the former question.

We can do the same, more or less, with a Fortran code and numpy F2PY. So we could start by defining subroutines phi_l and phi_t whose definitions are stored in file phi_for.f90:

subroutine phi_l(v, v_min, v_max, rate)
  real(kind=8), intent(in) :: v, v_min, v_max
  real(kind=8), intent(out) :: rate
  if (v < v_min) then
     rate = 0.0_8
  else if (v > v_max) then
     rate = 1.0_8
  else
     rate = (v-v_min)/(v_max-v_min)
  end if
end subroutine phi_l

subroutine phi_t(v, offset, slope, rate)
  real(kind=8), intent(in) :: v, offset, slope
  real(kind=8), intent(out) :: rate
  rate = (1.0_8 + tanh((v-offset)*slope))*0.5_8   
end subroutine phi_t

We compile it to make it usable in Python with:

python3 -m numpy.f2py -c phi_for.f90 -m phi_for

We test it against our initial version in the "linear" case:

import phi_for
all([phi_for.phi_l(i/25,0,1) == phi_linear(i/25,0,1) for i in range(26)])

True

And in the hyperbolic tangent one:

all([phi_for.phi_t(i/5,0,1) == phi_tanh(i/5,0,1) for i in range(26)])

True

Make a time estimation for phi_for.phi_l and phi_for.phi_t like in the former questions.

• The absolute value of the weight should be coded by the arrows thickness (see penwidth in the doc.).
• The value of the delay should be coded by the length of the edge (see weight in the doc.).
• The sign of the weight should be coded by the arrowhead: a positive sign leading to "–>" and a negative one to "–<" (see arrowType in the doc.).

If we have a three neurons, N0, N1, N2, graph such that:

• N0 inhibits N1 with weight -2 and delay 1 and N2 with weight -1 and delay 3.
• N1 excites N0 with weight 1 and delay 3 and N2 with weight 2 and delay 2.
• N2 excites N1 with weight 3 and delay 1.

The content of the string should be someting like:

digraph {
N1 -> N0 [arrowhead=normal penwidth=1 weight=0];
N0 -> N1 [arrowhead=inv penwidth=2 weight=7];
N2 -> N1 [arrowhead=normal penwidth=3 weight=3];
N0 -> N2 [arrowhead=inv penwidth=1 weight=0];
N1 -> N2 [arrowhead=normal penwidth=2 weight=3];
}

Notice that the representation of the weights magnitude (edge thickness) and delays (edge length) is only qualitative.

Create a Network class that combines compatible objects of classes SynapticGraph and NetworkDynamics. Here "compatible" means that they should have the same value for their _nb_types / _nb_kinds attributes. The class instances should have everything required to run the simulation. In particular the current state of the system should be managed by the instance and its methods. What makes the state will depend on the network structure—what is contained in the SynapticGraph instance—and on its dynamics—what is contained in the NetworkDynamics instance—. For instance, the state space of simplest network introduced in Sec. 1.4 consist in \(\mathbb{R}^n\) (\(n\), introduced at the beginning of Sec. 1.1, is here the value of attribute _size of the SynapticGraph instance) and its elements were defined by Eq. \eqref{eq:bigV}:

The dynamics is specified by Eq. \eqref{eq:simple-discrete-V-evolution}.

If a non trivial refractory period is introduced (Sec. 1.5.1), the state space must be extended to \((\mathbb{R}^n \times \mathbb{Z}^n)\). We must consider the joint process: \(\left(\mathbb{V}_t,\mathbb{L}_t\right)_{t \geq 0 }\), where the transitions of \(\mathbb{V}_t\) are given by Eq. \eqref{eq:discrete-V-evolution-with-rp} and the transition of \(\mathbb{L}_t\) are given by Eq. \eqref{eq:discrete-L-evolution}. \(\mathbb{L}_t\) is a vector valued random variable defined in Eq. \eqref{eq:bigL}:

If synaptic delays or non geometric, but compact, leakages/couplings are considered (Sec. 1.5.2 and 1.6.2) the state space must be extended further to: \(\mathbb{R}^n \times \mathbb{Z}^n \times \{0,1\}^{n \times (\tau_{max}+t_{max}+1)}\), and we must consider the joint process \((\mathbb{S}_t)_{t \ge 0} \equiv \left(\mathbb{V}_t,\mathbb{L}_t, \mathbb{X}_t\right)_{t \geq 0 }\), where \(\mathbb{X}_t\) is defined by Eq. \eqref{eq:bigX} or \eqref{eq:bigbigX}:

Your Network objects should therefore have the following extra "private" attribute (they will also have the attributes of the two classes they are derived from):

_t_now

the time (step) reached by the simulation.

_past_len

the length of the past to consider, SynapticGraph._tau_max + NetworkDynamics._t_max

_V

an array (of module array) of doubles with _size elements, it will contain \(\mathbb{V}_{t}\).

_L

an array of signed long integers with either no elements (if _r_max is 1) or with _size elements, it will contain \(\mathbb{L}_{t}\).

_X

a tuple of unsigned char arrays; the tuple will be of null length if there is no synaptic delay and no leakage or of length _size. Its elements when the tuple has an actual content will be arrays of size _t_max + _tau_max + 1. This attribute will contain \(\mathbb{X}_{t}\).

This is simply done with:

import timeit
left = 0
right = 1
print('phi_linear call:')
print(timeit.timeit('[phi_linear(i/25,left,right) for i in range(25)]',globals=globals()))
print('phi_linear2 call:')
print(timeit.timeit('[phi_linear2(i/25,left,right) for i in range(25)]',globals=globals()))
print('Direct estimation:')
print(timeit.timeit('[min(1,max(0,(i/25-left)/(right-left))) for i in range(25)]',globals=globals()))

phi_linear call:
6.551414999004919
phi_linear2 call:
10.65268645301694
Direct estimation:
8.990929204970598

The locally linear functions are simply compared with:

import timeit
left = 0
right = 1
print('phi_linear call:')
print(timeit.timeit('[phi_linear(i/25,left,right) for i in range(25)]',globals=globals()))
print('_phi_clib.phi_l call:')
print(timeit.timeit('[_phi_clib.phi_l(i/25,left,right) for i in range(25)]',globals=globals()))

phi_linear call:
6.389747657987755
_phi_clib.phi_l call:
22.24291676498251

The hyperbolic tangent based functions give:

import timeit
left = 0
right = 1
print('phi_tanh call:')
print(timeit.timeit('[phi_tanh(i/5,left,right) for i in range(25)]',globals=globals()))
print('_phi_clib.phi_t call:')
print(timeit.timeit('[_phi_clib.phi_t(i/5,left,right) for i in range(25)]',globals=globals()))

phi_tanh call:
30.951816001033876
_phi_clib.phi_t call:
24.338442599982955

This is simply done with:

import timeit
left = 0
right = 1
print('phi_for.phi_l call:')
print(timeit.timeit('[phi_for.phi_l(i/25,left,right) for i in range(25)]',globals=globals()))
print('phi_for.phi_t call:')
print(timeit.timeit('[phi_for.phi_t(i/5,left,right) for i in range(25)]',globals=globals()))

phi_for.phi_l call:
7.309892535035033
phi_for.phi_t call:
8.534695383976214

The class skeleton is:

class Neuron:
    <<Neuron-doc>>
    __slots__ = '_idx', '_kind', '_pre'
    <<Neuron-init>>
    <<Neuron-len>>
    <<Neuron-getitem>>
    <<Neuron-vset>>
    @property
    def idx(self):
        return self._idx
    @property
    def kind(self):
        return self._kind

The class docstring:

"""A neuron representing class.

A Neuron object has the following 'private' attributes:
- _idx: its index (positive or null integer) that uniquely identifies it.
- _kind: its kind/type (positive or null integer)
- _pre: a tuple with as many elements as the neuron receives presynaptic inputs.
        Each element is itself a tuple with 4 elements:
        - pre_idx: the index of the presynaptic neuron
        - pre_weight: the weight of the synapse (a real number)
        - pre_delay: the synaptic delay (positive or null integer)
        - pre_kind: the presynaptic neuron kind/type.

The class has the following methods:
- __init__: the constructor (not supposed to be called directly by the user).
- idx: returns the Neuron's _idx attribute (can be called with obj.idx).  
- kind: returns the Neuron's _kind attribute (can be called with obj.kind).
- len: returns the number of presynaptic neurons.
- []: Neuron[idx] returns the tuple of presynaptic neuron 'idx'
- vset: returns the list of indices of presynaptic neurons
        (potentially depending on their kind/type).
"""

We define here the special method __init__ that gets called upon class instance creation. We give it as parameters idx (default 0), kind (default 0) and pre (default (), an empty tuple):

def __init__(self,idx=0, kind=0, pre=()):
    """Initializes a Neuron.

    Parameters
    ----------
    idx: positive or null integer, the neuron index within the network.
    kind: positive or null integer, the neuron kind/type.
    pre: an iterable of tuples (pre_idx, pre_weight, pre_delay, pre_kind)
         where 'pre_idx' is the index of the presynaptic neuron,
               'pre_weight' is the weight of the synapse,
               'pre_delay', is the synaptic delay,
               'pre_kind', is the presynaptic kind/type.
         All the 'pre_idx' should be different.

    Returns
    -------
    Nothing but attributes: _idx, _kind and _pre are initialized.
    """
    from collections.abc import Iterable
    if not isinstance(idx,int):
        raise TypeError('idx must be an integer.')
    if idx < 0:
        raise ValueError('idx must be >= 0.')
    self._idx = idx
    if not isinstance(kind,int):
        raise TypeError('kind must be an integer.')
    if kind < 0:
        raise ValueError('kind must be >= 0.')
    self._kind = kind
    if not isinstance(pre,Iterable):
        raise TypeError('pre must be an Iterable.')
    if len(pre) > 0:
        if not all([isinstance(p,tuple) for p in pre]):
            raise ValueError('pre must contain tuples.')
        if not all([len(p)==4 for p in pre]):
            raise ValueError('pre must contain tuples of length 4.')
        if not all([all([isinstance(p[0],int),isinstance(p[1],float) or isinstance(p[1],int),
                        isinstance(p[2],int),isinstance(p[3],int)]) for p in pre]):
            raise ValueError('pre must contain (int,float,int,int) tuples.')
        if any([any([p[0] < 0, p[2] < 0, p[3] < 0]) for p in pre]):
            raise ValueError('pre must contain tuples with positive or null first, third and fourth elements.')
        if len(set([p[0] for p in pre])) != len(pre):
            raise ValueError('All first elements of the pre tuples should be different.')
    self._pre = tuple([t for t in pre])

This block defines the special method __len__ that makes the syntax len(obj) meaningful for an object of class Neuron by returning its number of presynaptic partners:

def __len__(self):
    """Returns the number of presynaptic partners."""
    return len(self._pre)

This block defines the special method __getitem__ that makes the syntax obj[i] meaningful for an object of class Neuron:

def __getitem__(self,idx):
    """Returns the tuple of presynaptic neuron 'idx'."""
    return self._pre[idx]

Naive instantiation:

n1 = Neuron()
len(n1)

0

More interesting a neuron with 3 presynaptic partners:

n2 = Neuron(1,0,[(0,0.5,1,0),(2,-0.5,2,1),(3,0.25,0,0)])
len(n2)

3

Let us look at the synaptic weights:

[t[1] for t in n2]

[0.5, -0.5, 0.25]

We can check vset with any of those:

n2.vset() # (0, 2, 3)
n2.vset(0) # (0, 3)
n2.vset(1) # (2,) 
n2.vset(2) # ()

Let us check that our checks work with any of:

n3 = Neuron(-1) # ValueError
n3 = Neuron(2,-1) # ValueError
n3 = Neuron(2,1,(0)) # TypeError
n3 = Neuron(2,1,(0,0)) # ValueError
n3 = Neuron(2,1,((0,),(0,))) # ValueError
n3 = Neuron(2,1,((0,1,2,3),(-1,1,2,3))) # ValueError
n3 = Neuron(1,0,[(0,0.5,1,0),(2,-0.5,2,1),(2,0.25,0,0)]) # ValueError

The class skeleton is:

class SynapticGraph:
    <<SynapticGraph-doc>>
    <<SynapticGraph-init>>
    <<SynapticGraph-getitem>>
    <<SynapticGraph-len>>
    <<SynapticGraph-dot>>
    @property
    def nb_types(self):
        return self._nb_types
    @property
    def tau_max(self):
        return self._tau_max

The class docstring:

"""A Synaptic graph class.

Objects of this class are essentially tuples of Neuron objects.
SynapticGraph objects have th dollowing 'private' attributes:
- _size: a postive integer, the nuber of neurons in the graph/network.
- _nb_types: a postive integer, the nuber of neuron kinds/types.
- _tau_max: a positive or null integer, the largest synaptic delay.
- _graph: a tuple of '_size' Neuron objects.

The following methods are defined:
- __init__: the class constructor (meant to be called by the user).
- nb_types: returns the value of attribute _nb_types (called with obj.nb_types).
- tau_max: returns the value of attribute _tau_max (called with obj.tau_max).  
- len: returns the number of neurons in the graph.
- []: SynapticGraph[idx] returns the corresponding Neuron object.
- dot: returns a string suitable for generating a graphical representation
       of the graph with GraphViz.
""" 

We define here the special method __init__ that gets called upon class instance creation. We give it as parameters weights and delays (optional).

def __init__(self, weights, delays=[], kinds=[], homogeneity=True):
    """Initializes instances of class SynapticGraph.

    Parameters
    ----------
    weights: a list of sub-lists. There are as many sub-lists as
             nodes/neurons in the directed graph. Each sub-list
             is as long as the parent list. The elements
             are real. The first index correspond to the pre-
             synaptic neuron, the second, to the postsynaptic.
    delays: an optional list of sub-lists. There are as many sub-
            lists as nodes/neurons in the directed graph. Each sub-
            list is as long as the parent list. The elements are
            positive or null integers. The first index correspond
            to the presynaptic neuron, the second, to the postsynaptic.
            If not specified all delays are 0.
    kinds: an optional list of postive or null integer specifying the
           kind/type of each neuron. There as many elements as neurons.
           If not specified, all the neurons are of the same kind.
    homogeneity: a Boolean, if True (default) the code checks that
                 all neurons of the same kind form synapses with the
                 same sign.

    Returns
    -------
    Nothing but attributes: size and graph get initialized.
    """
    from collections.abc import Iterable
    import numbers
    # Check weights
    <<SynapticGraph-check-weights>> 
    # Check delays
    <<SynapticGraph-check-delays>>
    <<SynapticGraph-check-kinds>>
    <<SynapticGraph-check-homogeneity>>
    <<SynapticGraph-graph-attribute>>
    

This block defines the special method __getitem__ that makes the syntax myobj[i] meaningful for an object of class SynapticGraph:

def __getitem__(self,idx):
    """Returns the tuple of presynaptic neuron triplets to neuron 'idx'."""
    return self._graph[idx]

This block defines the special method __len__ that makes the syntaxe len(myobj) meaningful for an object of class SynapticGraph by returning the network size:

def __len__(self):
    """Returns the network size."""
    return self._size

This block defines the dot method:

def dot(self,labels=[]):
    """Returns a string corresponding to the content of a DOT
    file suitable for generating a graphical representation.

    Parameters
    -----------
    labels: a list of strings, each element will be the label of
            one node/neuron of the graph. There should be as
            many elements as there are nodes/neurons. If left
            unspecified the labels are N0, N1, ...

    Returns
    --------
    A string.
    """
    if labels == []:
        labels = ["N"+str(i) for i in range(len(self))]
    txt = "digraph {\n"
    for post in range(len(self)):
        neuron = self[post]
        for pre_idx in range(len(neuron)):
            idx = neuron[pre_idx][0]
            line = "N"+str(idx)+" -> N"+str(post)+" ["
            if neuron[pre_idx][1] > 0:
                line += "arrowhead=normal "
            else:
                line += "arrowhead=inv "
            line += "penwidth="+str(abs(round(neuron[pre_idx][1])))+" "
            if self._tau_max > 0:
                w = 10*(1-neuron[pre_idx][2]/self.tau_max)
                line += "weight="+str(round(w))
            line += "];\n"
            txt += line
    txt += "}\n"
    return txt

Using the simple 3 neurons graph of the above example we do:

sg3 = SynapticGraph([[0,-2,-1],[1,0,2],[0,3,0]],
                    [[0,1,3],[3,0,2],[0,1,0]],
                    [1,0,0])
[tpl for tpl in sg3[1]]

[(0, -2, 1, 1), (2, 3, 1, 0)]

We can reproduce the previous graphical representation with:

with open("test-dot-for-SynapticGraph.dot","w") as f:
    f.write(sg3.dot())

Followed by the following instructions from the command line assuming that dot has been installed:

dot -Tpng test-dot-for-SynapticGraph.dot > test-dot-for-SynapticGraph.png

The class skeleton is:

class NeuronalDynamics:
    <<NeuronalDynamics-doc>>
    <<phi_linear>>
    <<phi_sigmoid>>
    <<NeuronalDynamics-init>>
    <<NeuronalDynamics-r>>
    <<NeuronalDynamics-phi>>
    <<NeuronalDynamics-g_self>>
    <<NeuronalDynamics-g_interaction>>
    @property
    def nb_kinds(self):
        return self._nb_kinds
    @property
    def r_max(self):
        return self._r_max
    @property
    def t_max(self):
        return self._t_max

The class docstring:

"""Required properties for specifying network dynamics.

The class 'private' attributes are:
- _nb_kinds: the number of neuron kinds/types (positive integer).
- _phi_tpl: a tuple of lenghth '_nb_kinds'. The spiking probability
            function of each neuronal kind/type.
- _rho_tpl: a tuple of lenghth 0 or '_nb_kinds'. In the latter case,
            each tuple element, rho, is such that 0 < rho < 1 and
            specifies a 'geometrical' leakage; that is, in the abscence
            of input the membrane potentials at two successive steps
            satisfy: v[t+1] = rho*v[t]
- _r_max: the longest refractory period (positive integer).
- _r_tpl: a tuple of length 0 (if '_r_max'=1) or  '_nb_kinds'. In the latter case,
          the refractory periods of each neuronal kind/type.
- _g_k: a tuple of tuples describing the self_interactions of the neurons of
        each neuronal kind/type. Can be of null length in the absence of
        self-interaction.
- _g_2k: a tuple of tuples of tuples describing the interactions between two
         neurons of two given kinds/types. Can be of null length if there is
         no leakage.

The class contains two functions that can be used to build spiking probability
functions:
- phi_linear
- phi_sigmoid

The class defines the following methods:
- r: NeuronalDynamics.r(idx) returns the refractory period of kind 'idx'.
- phi: NeuronalDynamics.phi(idx) returns the spiking probability function
       of kind 'idx'.
- g_self: NeuronalDynamics.g_self(idx,t) returns the self-effect of a neuron 
          of kind 'idx', 't' steps after it spiked.
- g_interaction: NeuronalDynamics.g_self(pre,post,t) returns the effect
                 of a spike from a neuron of kind 'pre' 't' steps after the
                 spike in a neuron of kind 'post'.
- nb_kinds: returns the value of attribute _nb_kinds (called with obj.nb_kinds).
- r_max: returns the value of attribute _r_max (called with obj.r_max).
- t_max: returns the value of attribute _t_max (called with obj.t_max).
"""

This methods returns the refractory period of kind/type idx:

def r(self,idx):
    """Refractory period of kind/type 'idx'."""
    if self.r_max == 1: return 1
    if not (0 <= idx < self.nb_kinds):
        raise ValueError('Must have 0 <= idx < ' + str(self.nb_kinds) + '.')
    else:
        return self._r_tpl[idx]
        

This method returns the spiking probability function of kind/type idx:

def phi(self,idx):
    """Spiking probability function of kind/type 'idx'."""
    if not (0 <= idx < self.nb_kinds):
        raise ValueError('Must have 0 <= idx < ' + str(self.nb_kinds) + '.')
    else:
        return self._phi_tpl[idx]
    

This method returns the "self-interaction" of kind idx at delay of t time steps after a spike:

def g_self(self,idx,t):
    """Self interaction of kind/type 'idx' at 't'."""
    if self.t_max == 0: return 0.0
    if not (0 <= idx < self.nb_kinds):
        raise ValueError('Must have 0 <= idx < ' + str(self.nb_kinds) + '.')
    else:
        if not 0 <= t <= self.t_max: return 0.0
        else:
            return self._g_k[idx][t]
        

This method returns the effect of a spike from a neuron of kind/type pre_idx in a neuron of kind/type post_idx after a delay of t time steps:

def g_interaction(self, pre_idx, post_idx, t):
    """Interaction of kind 'pre_idx' with 'post_idx' at 't'."""
    if self.t_max == 0: return 1.0
    if post_idx == pre_idx: return self.g_self(pre_idx,t)
    if not (0 <= pre_idx < self.nb_kinds):
        raise ValueError('Must have 0 <= pre_idx < ' + str(self.nb_kinds) + '.')
    if not (0 <= post_idx < self.nb_kinds):
        raise ValueError('Must have 0 <= post_idx < ' + str(self.nb_kinds) + '.')
    if not 0 <= t <= self.t_max: return 0.0
    else:
        return self._g_2k[pre_idx][post_idx][t]
 

The class skeleton is:

class Network(SynapticGraph,NeuronalDynamics):
    <<Network-doc>>
    array = __import__('array')
    <<Network-init>>
    @property
    def t_now(self):
        return self._t_now
    @property
    def past_len(self):
        return self._past_len
    <<Network-spike-prob>>
    <<Network-update_nl_simple>>

Class docstring;

def __init__(self, SG_obj, ND_obj, V0 = None, L0 = None, X0 = None):
    <<Network-init-doc>>
    from collections.abc import Iterable
    <<Network-check-SG_obj-ND_obj>>
    <<Network-set-inherited-attributes>>
    self._past_len = ND_obj.t_max + SG_obj.tau_max
    <<Network-set-_V-attribute>>
    <<Network-set-_L-attribute>>
    <<Network-set-_X-attribute>>
    self._t_now = 0

Method spike_prob returns the spiking probability of each neuron of the network as a list. The property decorator makes it callable with obj.spike_prob.

@property
def spike_prob(self):
    """Spiking probability of each neuron of the Network."""
    return [self.phi(self[i].kind)(self._V[i]) for i in range(len(self))]

This block defines the simplest update method for the non-leakage case without refractory period and without synaptic delay corresponding to Eq. \eqref{eq:simple-discrete-V-evolution}:

def update_nl_simple(self,spike_or_not):
    """Update Network state depending on Boolean vector 'spike_or_not'
    for the simplest model: no leakge, no refractory period, no delay."""
    self._t_now += 1
    #breakpoint()
    for n_idx in range(len(self)):
        if spike_or_not[n_idx]: # neuron n_idx spiked
            self._V[n_idx] = 0
        else: # neuron n_idx did not spike
            for pre in self[n_idx]:
                if spike_or_not[pre[0]]: # neuron pre[0] spiked
                    self._V[n_idx] += pre[1]