Source code:
https://bitbucket.org/annarchy/annarchy
Documentation:
https://annarchy.readthedocs.io/en/stable/
Forum:
https://groups.google.com/forum/#!forum/annarchy
Notebooks used in this tutorial:
Installation guide: https://annarchy.readthedocs.io/en/stable/intro/Installation.html
From pip:
pip install ANNarchy
From source:
git clone https://bitbucket.org/annarchy/annarchy.git
cd annarchy
python setup.py install
Requirements (Linux and MacOS):
Simulation of both rate-coded and spiking neural networks.
Only local biologically realistic mechanisms are possible (no backpropagation).
Equation-oriented description of neural/synaptic dynamics (à la Brian).
Code generation in C++, parallelized using OpenMP on CPU and CUDA on GPU (MPI is coming).
Synaptic, intrinsic and structural plasticity mechanisms.
from ANNarchy import *
setup(dt=1.0)
neuron = Neuron(...) # Create neuron types
stdp = Synapse(...) # Create synapse types for transmission and/or plasticity
pop = Population(1000, neuron) # Create populations of neurons
proj = Projection(pop, pop, 'exc', stdp) # Connect the populations
proj.connect_fixed_probability(weights=0.0, probability=0.1)
compile() # Generate and compile the code
m = Monitor(pop, ['spike']) # Record spikes
simulate(1000.) # Simulate for 1 second
data = m.get('spike') # Retrieve the data and plot it
ESN rate-coded neurons typically follow first-order ODEs:
$$ \tau \frac{dx(t)}{dt} + x(t) = \sum w^\text{in} , r^\text{in}(t) + g , \sum w^\text{rec} , r(t) + \xi(t) $$
$$ r(t) = \tanh(x(t)) $$
from ANNarchy import *
ESN_Neuron = Neuron(
parameters = """
tau = 30.0 # Time constant
g = 1.0 : population # Scaling
noise = 0.01 : population # Noise amplitude
""",
equations="""
tau * dx/dt + x = sum(in) + g * sum(exc) + noise * Uniform(-1, 1) : init=0.0
r = tanh(x)
"""
)
parameters = """
tau = 30.0 # Time constant
g = 1.0 : population # Scaling
noise = 0.01 : population # Noise amplitude
"""
All parameters used in the equations must be declared in the Neuron definition.
Parameters can have one value per neuron in the population (default) or be common to all neurons (flag population
or projection
).
Parameters and variables are double floats by default, but the type can be specified (int
, bool
).
equations="""
tau * dx/dt + x = sum(in) + g * sum(exc) + noise * Uniform(-1, 1) : init=0.0
r = tanh(x)
"""
Variables are evaluated at each time step in the order of their declaration, except for coupled ODEs.
Variables can be updated with assignments (=
, +=
, etc) or by defining first order ODEs.
The math C library symbols can be used (tanh
, cos
, exp
, etc).
Initial values at $t=0$ can be specified with init
(default: 0.0).
Lower/higher bounds on the values of the variables can be set with the min
/max
flags:
r = x : min=0.0 # ReLU
Additive noise can be drawn from several distributions, including Uniform
, Normal
, LogNormal
, Exponential
, Gamma
…
The output variable of a rate-coded neuron must be r
.
First-order ODEs are parsed and manipulated using sympy
:
# All equivalent:
tau * dx/dt + x = 0.0
tau * dx/dt = - x
dx/dt = (-x)/tau
Several numerical methods are available (https://annarchy.readthedocs.io/en/stable/manual/NumericalMethods.html):
Explicit (forward) Euler (default): tau * dx/dt + x = 0.0 : init=0.0, explicit
Implicit (backward) Euler: tau * dx/dt + x = 0.0 : init=0.0, implicit
Exponential Euler (exact for linear ODE): tau * dx/dt + x = 0.0 : init=0.0, exponential
Midpoint (RK2): tau * dx/dt + x = 0.0 : init=0.0, midpoint
Event-driven (spiking synapses): tau * dx/dt + x = 0.0 : init=0.0, event-driven
ODEs are solved concurrently, instead of sequentially for assignments:
# I is updated
I = sum(exc) - sum(inh) + b
# u and v are solved concurrently using the current of I
tau * dv/dt + v = I - u
tau * du/dt + u = v
# r uses the updated value of v
r = tanh(v)
The order of the equations therefore matters a lot.
A single variable can only be updated once in the equations
field.
Populations are creating by specifying a number of neurons and a neuron type:
pop = Population(1000, ESN_Neuron)
For visualization purposes or when using convolutional layers, a tuple geometry can be passed instead of the size:
pop = Population((100, 100), ESN_Neuron)
All parameters and variables become attributes of the population (read and write) as numpy arrays:
pop.tau = np.linspace(20.0, 40.0, 1000)
pop.r = np.tanh(pop.v)
Single neurons can be individually modified, if the population
flag was not set:
pop[10].r = 1.0
Slices of populations are called PopulationView
and can be addressed separately:
pop = Population(1000, ESN_Neuron)
E = pop[:800]
I = pop[800:]
Projections link two populations (or views) in a uni-directional way.
proj_exc = Projection(E, pop, 'exc')
proj_inh = Projection(I, pop, 'inh')
Each target ('exc', 'inh', 'AMPA', 'NMDA', 'GABA'
) can be defined as needed and will be treated differently by the post-synaptic neurons.
The weighted sum of inputs for a specific target is accessed in the equations by sum(target)
:
equations="""
tau * dx/dt + x = sum(exc) - sum(inh)
r = tanh(x)
"""
It is therefore possible to model modulatory effects, divisive inhibition, etc.
Projections must be populated with a connectivity matrix (who is connected to who), a weight w
and optionally a delay d
(uniform or variable).
Several patterns are predefined:
proj.connect_all_to_all(weights=Normal(0.0, 1.0), delays=2.0, allow_self_connections=False)
proj.connect_one_to_one(weights=1.0, delays=Uniform(1.0, 10.0))
proj.connect_fixed_number_pre(number=20, weights=1.0)
proj.connect_fixed_number_post(number=20, weights=1.0)
proj.connect_fixed_probability(probability=0.2, weights=1.0)
proj.connect_gaussian(amp=1.0, sigma=0.2, limit=0.001)
proj.connect_dog(amp_pos=1.0, sigma_pos=0.2, amp_neg=0.3, sigma_neg=0.7, limit=0.001)
But you can also load Numpy arrays or Scipy sparse matrices. Example for synfire chains:
w = np.array([[None]*pre.size]*post.size)
for i in range(post.size):
w[i, (i-1)%pre.size] = 1.0
proj.connect_from_matrix(w)
w = lil_matrix((pre.size, post.size))
for i in range(pre.size):
w[pre.size, (i+1)%post.size] = 1.0
proj.connect_from_sparse(w)
Once all populations and projections are created, you have to generate to the C++ code and compile it:
compile()
You can now manipulate all parameters/variables from Python thanks to the Cython bindings.
A simulation is simply run for a fixed duration with:
simulate(1000.) # 1 second
You can also run a simulation until a criteria is filled, check:
https://annarchy.readthedocs.io/en/stable/manual/Simulation.html#early-stopping
By default, a simulation is run in C++ without interaction with Python.
You may want to record some variables (neural or synaptic) during the simulation with a Monitor
:
m = Monitor(pop, ['v', 'r'])
n = Monitor(proj, ['w'])
After the simulation, you can retrieve the recordings with:
recorded_v = m.get('v')
recorded_r = m.get('r')
recorded_w = n.get('w')
Warning: calling get()
flushes the array.
Warning: recording projections can quickly fill up the RAM (see Dendrites).
Link to the Jupyter notebook on github: RC.ipynb
Spiking neurons must also define two additional fields:
spike
: condition for emitting a spike.
reset
: what happens after a spike is emitted (at the start of the refractory period).
A refractory period in ms can also be specified.
Example of the Leaky Integrate-and-Fire:
LIF = Neuron(
parameters="""
tau = 20.
E_L = -70.
v_T = 0.
v_r = -58.
I = 50.0
""",
equations="""
tau * dv/dt = (E_L - v) + I : init=E_L
""",
spike=" v >= v_T ",
reset=" v = v_r ",
refractory = 2.0
)
A pre-synaptic spike arriving to a spiking neuron increase the conductance g_target
(e.g. g_exc
or g_inh
, depending on the projection).
LIF = Neuron(
parameters="...",
equations="""
tau * dv/dt = (E_L - v) + g_exc - g_inh
""",
spike=" v >= v_T ",
reset=" v = v_r ",
refractory = 2.0
)
Each spike increments g_target
from the synaptic efficiency w
of the corresponding synapse.
g_target += w
This defines an instantaneous model of synaptic transmission.
For exponentially-decreasing or alpha-shaped synapses, ODEs have to be introduced for the conductance/current.
The exponential numerical method should be preferred, as integration is exact.
LIF = Neuron(
parameters="...",
equations="""
tau * dv/dt = (E_L - v) + g_exc + alpha_exc # exponential or alpha
tau_exc * dg_exc/dt = - g_exc : exponential
tau_exc * dalpha_exc/dt = exp((tau_exc - dt/2.0)/tau_exc) * g_exc
- alpha_exc : exponential
""",
spike=" v >= v_T ",
reset=" v = v_r ",
refractory = 2.0
)
Link to the Jupyter notebook on github: AdEx.ipynb
$$ \tau , \frac{dv}{dt} = (E_L - v) + \delta_T , \exp \frac{v-v_T}{\delta_T} + I - w $$ $$ \tau_w , \frac{dw}{dt} = a , (v - E_L) - w $$
AdEx = Neuron(
parameters="""
tau = 20.
E_L = -70.
v_T = -50. ; v_r = -58.
delta_T = 2.0
a = 0.2 ; b = 0.
tau_w = 30.
I = 50.0
""",
equations="""
tau * dv/dt = (E_L - v) + delta_T * exp((v-v_T)/delta_T) + I - w : init=E_L
tau_w * dw/dt = a * (v - E_L) - w : init=0.0
""",
spike=" v >= 0.0 ",
reset=" v = v_r ; w += b ",
refractory = 2.0
)
Synapses can also implement equations that will be evaluated after each neural update.
IBCM = Synapse(
parameters = """
eta = 0.01 : projection
tau = 2000.0 : projection
""",
equations = """
tau * dtheta/dt + theta = post.r^2 : postsynaptic, exponential
dw/dt = eta * post.r * (post.r - theta) * pre.r : min=0.0, explicit
""",
psp = " w * pre.r"
)
The synaptic efficiency (weight) must be w
.
Each synapse can access pre- and post-synaptic variables with pre.
and post.
.
The postsynaptic
flag allows to do computations only once per post-synaptic neurons.
psp
optionally defines what will be summed by the post-synaptic neuron (e.g. psp = "w * log(pre.r)"
).
The synapse type just has to be passed to the Projection:
proj = Projection(inp, pop, 'exc', IBCM)
Synaptic variables can be accessed as lists of lists for the whole projection:
proj.w
proj.theta
or for a single post-synaptic neuron (Dendrite
):
proj[10].w
Link to the Jupyter notebook on github: Miconi.ipynb
Spiking synapses can define a pre_spike
field, defining what happens when a pre-synaptic spike arrives at the synapse.
g_target
is an alias for the corresponding post-synaptic conductance: it will be replaced by g_exc
or g_inh
depending on how the synapse is used.
By default, a pre-synaptic spike increments the post-synaptic conductance from w
: g_target += w
STP = Synapse(
parameters = """
tau_rec = 100.0 : projection
tau_facil = 0.01 : projection
U = 0.5
""",
equations = """
dx/dt = (1 - x)/tau_rec : init = 1.0, event-driven
du/dt = (U - u)/tau_facil : init = 0.5, event-driven
""",
pre_spike="""
g_target += w * u * x
x *= (1 - u)
u += U * (1 - u)
"""
)
post_spike
similarly defines what happens when a post-synaptic spike is emitted.
STDP = Synapse(
parameters = """
tau_plus = 20.0 : projection ; tau_minus = 20.0 : projection
A_plus = 0.01 : projection ; A_minus = 0.01 : projection
w_min = 0.0 : projection ; w_max = 1.0 : projection
""",
equations = """
tau_plus * dx/dt = -x : event-driven # pre-synaptic trace
tau_minus * dy/dt = -y : event-driven # post-synaptic trace
""",
pre_spike="""
g_target += w
x += A_plus * w_max
w = clip(w + y, w_min , w_max)
""",
post_spike="""
y -= A_minus * w_max
w = clip(w + x, w_min , w_max)
""")
Standard populations (SpikeSourceArray
, TimedArray
, PoissonPopulation
, HomogeneousCorrelatedSpikeTrains
), OpenCV bindings.
Standard neurons:
Standard synapses:
Parallel simulations with parallel_run
.
Convolutional and pooling layers.
Hybrid rate-coded / spiking networks.
Structural plasticity.