Copyright (C) 2014 by Eka A. Kurniawan

eka.a.kurniawan(ta)gmail(tod)com

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; version 2 of the License.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

TODOs

  • Parallel Python-OpenCL Implementation: Implement spiking input neuron with synaptic delay using integer event.
  • Parallel Python-OpenCL Implementation: Implement spiking input neuron with synaptic delay using binary event.
  • Execution Time: Collect and plot execution time based on different numbers of outgoing synapses.
  • Memory Consumption: Calculate memory consumption for different outgoing synapse sizes using both integer and binary event.

Done:

  • Sequential Python Implementation: Implement spiking input neuron with synaptic delay using binary event.
  • Sequential Python Implementation: Implement spiking input neuron with synaptic delay using integer event.
  • Spiking Input Neuron Model with Synaptic Delay: Add spiking input neuron class.

Tested On

MacBook Pro Retina, Mid 2012 with OS X 10.9.1 (Mavericks)

Python version:


In [1]:
import sys
print("Python %d.%d.%d" % (sys.version_info.major, \
                           sys.version_info.minor, \
                           sys.version_info.micro))


Python 3.3.3

NumPy version:


In [2]:
import numpy as np
print("NumPy %s" % np.__version__)


NumPy 1.8.0

matplotlib version:


In [3]:
import matplotlib
import matplotlib.pyplot as plt
print("matplotlib %s" % matplotlib.__version__)


matplotlib 1.3.1

Execution and Settings

To run this IPython Notebook; open a console, go to notebook directory and execute following command.

ipython3-3.3 notebook --pylab inline

Settings required:


In [4]:
# Display graph in 'retina' format for Mac with retina display. Others, use PNG or SVG format.
%config InlineBackend.figure_format = 'retina'
#%config InlineBackend.figure_format = 'PNG'
#%config InlineBackend.figure_format = 'SVG'

Other imports:


In [5]:
import time
from random import randint
import os.path

Housekeeping Functions

Functions to save/load array to/from file.


In [6]:
def save_list(file_name, list_name):
    with open(file_name, 'w') as fh:
        fh.write("\n".join(list_name))
    
def load_list(file_name):
    with open(file_name, 'r') as fh:
        return [line.strip() for line in fh.readlines()]

Function to print event buffer with the current index location (marked with >>>).


In [7]:
# Print event in array of integers
def print_event_ie(cur_event_idx = None, events = [[]]):
    for event_id, event in enumerate(events):
        if event_id == cur_event_idx:
            print(">>> %s" % event)
        else:
            print("    %s" % event)
            
# Print event in array of binaries
def print_event_be(events = []):
    for event in events:
        print('{0:011b}'.format(event))

Spiking Input Neuron Model with Synaptic Delay

The implementation is taken from CSIM: A neural Circuit SIMulator website and modified accordingly.[1]


In [8]:
# incoming_synapses: incoming synapses ID
# outgoing_synapses: outgoing synapses ID
class Neuron:
    def __init__(self,
                 incoming_synapses = [],
                 outgoing_synapses = []):
        self.incoming_synapses = incoming_synapses
        self.outgoing_synapses = outgoing_synapses

# Spiking input neuron class with synaptic delay using integer event (ie).
# Spiking input neuron has no incoming synapses as the channel attached 
# to the neuron provides all input spikes.
# outgoing_delays: array of outgoing synaptic delays connected to 
#                  this particular neuron in time step
# channel: None ternimated array of simulation time when spikes happen
class SpikingInputNeuron_SynapticDelay_ie(Neuron):
    def __init__(self,
                 outgoing_synapses = [],
                 outgoing_delays = [],
                 channel = []):
        Neuron.__init__(self, [], outgoing_synapses)
        self.outgoing_delays = outgoing_delays
        self.channel = channel
        self.nextSpikeTime = channel[0]
        # Current channel index
        self.cur_idx = 0
        self.hasFired = 0

    # Get next value from channel
    def nextValue(self, t):
        self.hasFired = 0
        
        if self.nextSpikeTime == None:
            return
        
        if self.nextSpikeTime < t:
            self.hasFired = 1
            self.cur_idx += 1
            self.nextSpikeTime = channel[self.cur_idx]
    
    # Handle synaptic delay and update event buffer
    def propagateSpike(self, cur_event_idx, ldelayList, outgoing_events):
        for outgoing_synapse in self.outgoing_synapses:
            delay = self.outgoing_delays[outgoing_synapse]
            new_event_idx = cur_event_idx + delay
            if new_event_idx >= ldelayList:
                new_event_idx -= ldelayList
            outgoing_events[new_event_idx][outgoing_synapse] = 1
    
    # Function that consists of all processes to be run every iteration
    def nextState(self, t, cur_event_idx, ldelayList, outgoing_events):
        self.nextValue(t)
        if self.hasFired:
            self.propagateSpike(cur_event_idx, ldelayList, outgoing_events)

# Spiking input neuron class with synaptic delay using binary event (be).
# Spiking input neuron has no incoming synapses as the channel attached 
# to the neuron provides all input spikes.
# outgoing_delays: array of outgoing synaptic delays connected to 
#                  this particular neuron in bit location of time step
# channel: None ternimated array of simulation time when spikes happen
class SpikingInputNeuron_SynapticDelay_be(Neuron):
    def __init__(self,
                 outgoing_synapses = [],
                 outgoing_delays = [],
                 channel = []):
        Neuron.__init__(self, [], outgoing_synapses)
        self.outgoing_delays = outgoing_delays
        self.channel = channel
        self.nextSpikeTime = channel[0]
        # Current channel index
        self.cur_idx = 0
        self.hasFired = 0

    # Get next value from channel
    def nextValue(self, t):
        self.hasFired = 0
        
        if self.nextSpikeTime == None:
            return
        
        if self.nextSpikeTime < t:
            self.hasFired = 1
            self.cur_idx += 1
            self.nextSpikeTime = channel[self.cur_idx]
    
    # Handle synaptic delay and update event buffer
    def propagateSpike(self, outgoing_events):
        for outgoing_synapse in self.outgoing_synapses:
            delay = self.outgoing_delays[outgoing_synapse]
            outgoing_events[outgoing_synapse] += delay
    
    # Function that consists of all processes to be run every iteration
    def nextState(self, t, outgoing_events):
        self.nextValue(t)
        if self.hasFired:
            self.propagateSpike(outgoing_events)

Test spiking input neuron using integer event with input spike at 0.801 millisecond.


In [9]:
dt = 2e-4
Tsim = 30e-4

ttl_synapse = 10
# Outgoing synapse IDs
outgoing_synapses = [1,3,7]
# Outgoing delay (in time step)
outgoing_delays = [0,2,0,10,0,0,0,5,0,0]
# Get maximum outgoing delay to alocate the most optimized event buffer
ldelayList = np.max(outgoing_delays) + 1
# 2D array of outgoing events connected to this particular neuron in integer
outgoing_events = [ttl_synapse * [0] for i in range(ldelayList)]
# Input channel
channel = [8.01e-4, None]

# Construct a spiking input neuron
n1 = SpikingInputNeuron_SynapticDelay_ie(outgoing_synapses, \
                                         outgoing_delays, \
                                         channel)

t = 0.0
nSteps = int((Tsim / dt) + 0.5)
cur_event_idx = 0
for step in range(nSteps):
    t += dt
    n1.nextState(t, cur_event_idx, ldelayList, outgoing_events)
    
    if int(channel[0] / dt) == step:
        print_event_ie(cur_event_idx, outgoing_events)
    
    cur_event_idx += 1
    if (cur_event_idx >= ldelayList):
        cur_event_idx = 0


    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
>>> [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Test spiking input neuron using integer event with input spike at 1.601 milliseconds.


In [10]:
dt = 2e-4
Tsim = 30e-4

ttl_synapse = 10
# Outgoing synapse IDs
outgoing_synapses = [1,3,7]
# Outgoing delay (in time step)
outgoing_delays = [0,2,0,10,0,0,0,5,0,0]
# Get maximum outgoing delay to alocate the most optimized event buffer
ldelayList = np.max(outgoing_delays) + 1
# 2D array of outgoing events connected to this particular neuron in integer
outgoing_events = [ttl_synapse * [0] for i in range(ldelayList)]
# Input channel
channel = [16.01e-4, None]

# Construct a spiking input neuron
n1 = SpikingInputNeuron_SynapticDelay_ie(outgoing_synapses, \
                                         outgoing_delays, \
                                         channel)

t = 0.0
nSteps = int((Tsim / dt) + 0.5)
cur_event_idx = 0
for step in range(nSteps):
    t += dt
    n1.nextState(t, cur_event_idx, ldelayList, outgoing_events)
    
    if int((channel[0] / dt) + 0.5) == step:
        print_event_ie(cur_event_idx, outgoing_events)
    
    cur_event_idx += 1
    if (cur_event_idx >= ldelayList):
        cur_event_idx = 0


    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
>>> [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]

Test spiking input neuron using integer event with input spike at 2.601 milliseconds.


In [11]:
dt = 2e-4
Tsim = 30e-4

ttl_synapse = 10
# Outgoing synapse IDs
outgoing_synapses = [1,3,7]
# Outgoing delay (in time step)
outgoing_delays = [0,2,0,10,0,0,0,5,0,0]
# Get maximum outgoing delay to alocate the most optimized event buffer
ldelayList = np.max(outgoing_delays) + 1
# 2D array of outgoing events connected to this particular neuron in integer
outgoing_events = [ttl_synapse * [0] for i in range(ldelayList)]
# Input channel
channel = [26.01e-4, None]

# Construct a spiking input neuron
n1 = SpikingInputNeuron_SynapticDelay_ie(outgoing_synapses, \
                                         outgoing_delays, \
                                         channel)

t = 0.0
nSteps = int((Tsim / dt) + 0.5)
cur_event_idx = 0
for step in range(nSteps):
    t += dt
    n1.nextState(t, cur_event_idx, ldelayList, outgoing_events)
    
    if int(channel[0] / dt) == step:
        print_event_ie(cur_event_idx, outgoing_events)
    
    cur_event_idx += 1
    if (cur_event_idx >= ldelayList):
        cur_event_idx = 0


    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
>>> [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Test spiking input neuron using binary event with input spike at 2.601 milliseconds.


In [12]:
dt = 2e-4
Tsim = 30e-4

ttl_synapse = 10
# Outgoing synapse IDs
outgoing_synapses = [1,3,7]
# Outgoing delay (in time step)
outgoing_delays_ie = [0,2,0,10,0,0,0,5,0,0]
# Outgoing delay (in binary location)
outgoing_delays_be = [2 ** delay for delay in outgoing_delays_ie]
# Get maximum outgoing delay is not allowed to be greather than 31
# for 32-bit integer data type
ldelayList = np.max(outgoing_delays_ie) + 1
if ldelayList > 31:
    print("WARNING: Maximum outgoing delay is %s time step(s)" % ldelayList)
    print("         It is not allowed to be greather than 31 for 32-bit integer data type")
# array of outgoing events connected to this particular neuron in 
# integer of binary
outgoing_events = np.array([0 for i in range(ttl_synapse)])
# Input channel
channel = [26.01e-4, None]

# Construct a spiking input neuron
n1 = SpikingInputNeuron_SynapticDelay_be(outgoing_synapses, \
                                         outgoing_delays_be, \
                                         channel)

t = 0.0
nSteps = int((Tsim / dt) + 0.5)
for step in range(nSteps):
    t += dt
    n1.nextState(t, outgoing_events)

    if int(channel[0] / dt) == step:
        print_event_be(outgoing_events)
    
    # Shift outgoing event one to right
    outgoing_events = outgoing_events >> 1


00000000000
00000000100
00000000000
10000000000
00000000000
00000000000
00000000000
00000100000
00000000000
00000000000

Sequential Python Implementation

Sequential Implementation using Integer Event (seq_ie Function)


In [13]:
def run_seq_ie(nSteps, ldelayList, \
               outgoing_synapses, outgoing_delays, outgoing_events, \
               channel, \
               plot = False):
    # Buffer allocation time is not applicable for sequential implementation
    bat = np.nan
    
    # Construct a spiking input neuron
    n1 = SpikingInputNeuron_SynapticDelay_ie(outgoing_synapses, \
                                             outgoing_delays, \
                                             channel)
    
    t = 0.0
    cur_event_idx = 0
    tic = time.time()
    for step in range(nSteps):
        t += dt
        n1.nextState(t, cur_event_idx, ldelayList, outgoing_events)
        cur_event_idx += 1
        if (cur_event_idx >= ldelayList):
            cur_event_idx = 0
    toc = time.time()
    # Run time
    rt = toc - tic
    # Total time
    tt = rt

    if plot:
        print_event_ie(None, outgoing_events)
    
    return bat, rt, tt

Test spiking input neuron with input spike at 2.601 milliseconds.


In [14]:
dt = 2e-4
Tsim = 0.3

ttl_synapse = 10
# Outgoing synapse IDs
outgoing_synapses = [1,3,7]
# Outgoing delay (in time step)
outgoing_delays = [0,2,0,10,0,0,0,5,0,0]
# Get maximum outgoing delay to alocate the most optimized event buffer
ldelayList = np.max(outgoing_delays) + 1
# 2D array of outgoing events connected to this particular neuron in integer
outgoing_events = [ttl_synapse * [0] for i in range(ldelayList)]
# Input channel
channel = [26.01e-4, None]

nSteps = int((Tsim / dt) + 0.5)
seq_t = run_seq_ie(nSteps, ldelayList, \
                   outgoing_synapses, outgoing_delays, outgoing_events, \
                   channel, \
                   plot = True)


    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Sequential Implementation using Binary Event (seq_be Function)


In [15]:
def run_seq_be(nSteps, \
               outgoing_synapses, outgoing_delays, outgoing_events, \
               channel, \
               plot = False):
    # Buffer allocation time is not applicable for sequential implementation
    bat = np.nan
    
    # Construct a spiking input neuron
    n1 = SpikingInputNeuron_SynapticDelay_be(outgoing_synapses, \
                                             outgoing_delays, \
                                             channel)
    
    t = 0.0
    tic = time.time()
    for step in range(nSteps):
        t += dt
        n1.nextState(t, outgoing_events)
        # Shift outgoing event one to right
        outgoing_events = outgoing_events >> 1
    toc = time.time()
    # Run time
    rt = toc - tic
    # Total time
    tt = rt

    if plot:
        print_event_be(outgoing_events)
    
    return bat, rt, tt

Test spiking input neuron with input spike at 2.601 milliseconds.


In [16]:
dt = 2e-4
Tsim = 28.01e-4

ttl_synapse = 10
# Outgoing synapse IDs
outgoing_synapses = [1,3,7]
# Outgoing delay (in time step)
outgoing_delays_ie = [0,2,0,10,0,0,0,5,0,0]
# Outgoing delay (in binary location)
outgoing_delays_be = [2 ** delay for delay in outgoing_delays_ie]
# Get maximum outgoing delay is not allowed to be greather than 31
# for 32-bit integer data type
ldelayList = np.max(outgoing_delays_ie) + 1
if ldelayList > 31:
    print("WARNING: Maximum outgoing delay is %s time step(s)" % ldelayList)
    print("         It is not allowed to be greather than 31 for 32-bit integer data type")
# array of outgoing events connected to this particular neuron in 
# integer of binary
outgoing_events = np.array([0 for i in range(ttl_synapse)])
# Input channel
channel = [26.01e-4, None]

nSteps = int((Tsim / dt) + 0.5)
seq_t = run_seq_be(nSteps, \
                   outgoing_synapses, outgoing_delays_be, outgoing_events, \
                   channel, \
                   plot = True)


00000000000
00000000010
00000000000
01000000000
00000000000
00000000000
00000000000
00000010000
00000000000
00000000000

Execution Time

Following code executes serial and all parallel implementations on different number of outgoing synapses (1, 10, 100, 1 thousand, 10 thousand, 100 thousand, 1 million) and different OpenCL devices supported by this machine. The execution time collected from each implementation consists of buffer allocation time (bat), run time (rt) and total time (tt) of the two.


In [17]:
dt = 2e-4
Tsim = 0.3

ttl_synapses = [1, 10, 100, 1000, 10000, 100000, 1000000]

# Outgoing synaptic delay range. Both start and end are included.
osd_start = 2
osd_end = 11

# Input channel
channel = [10e-3, 15e-3, 20e-3, 25e-3, 30e-3, 35e-3, 40e-3, 45e-3, 50e-3, 55e-3] + \
          [100e-3] + \
          [120e-3, 130e-3, 140e-3, 150e-3, 155e-3, 160e-3, 165e-3, 170e-3] + \
          [250e-3, 255e-3, 260e-3, 265e-3, 270e-3, 280e-3, 290e-3] + \
          [None]

nSteps = int((Tsim / dt) + 0.5)
for ttl_synapse in ttl_synapses:
    print(ttl_synapse) #bar
    # Outgoing synapse IDs
    outgoing_synapses = [sid for sid in range(ttl_synapse)]
    # Outgoing delay (in time step)
    outgoing_delays_ie = [randint(osd_start, osd_end) for sid in range(ttl_synapse)]
    # Outgoing delay (in binary location)
    outgoing_delays_be = [2 ** delay for delay in outgoing_delays_ie]
    
    # Get maximum outgoing delay to alocate the most optimized event buffer
    ldelayList = np.max(outgoing_delays) + 1
    if ldelayList > 31:
        print("WARNING: Maximum outgoing delay is %s time step(s)" % ldelayList)
        print("         It is not allowed to be greather than 31 for 32-bit integer data type")
    # Alocate event buffer
    outgoing_events_ie = [ttl_synapse * [0] for i in range(ldelayList)]
    outgoing_events_be = np.array([0 for i in range(ttl_synapse)])
    
    seq_ie_t = run_seq_ie(nSteps, ldelayList, \
                          outgoing_synapses, outgoing_delays_ie, outgoing_events_ie, \
                          channel)
    print(seq_ie_t) #bar
    
    seq_be_t = run_seq_be(nSteps, \
                          outgoing_synapses, outgoing_delays_be, outgoing_events_be, \
                          channel)
    print(seq_be_t) #bar


1
(nan, 0.0031969547271728516, 0.0031969547271728516)
(nan, 0.008183002471923828, 0.008183002471923828)
10
(nan, 0.0036079883575439453, 0.0036079883575439453)
(nan, 0.009937047958374023, 0.009937047958374023)
100
(nan, 0.00824117660522461, 0.00824117660522461)
(nan, 0.012319087982177734, 0.012319087982177734)
1000
(nan, 0.0534818172454834, 0.0534818172454834)
(nan, 0.055780887603759766, 0.055780887603759766)
10000
(nan, 0.49008798599243164, 0.49008798599243164)
(nan, 0.4440939426422119, 0.4440939426422119)
100000
(nan, 4.807276964187622, 4.807276964187622)
(nan, 4.334025144577026, 4.334025144577026)
1000000
(nan, 47.88603591918945, 47.88603591918945)
(nan, 43.684969902038574, 43.684969902038574)

References