arbor_test_sde.py

# git clone --depth 1 --recursive --branch=SDE 'https://github.com/boeschf/arbor'
# cd arbor ; mkdir build ; cd build
# cmake -DCMAKE_INSTALL_PREFIX=$(realpath prefix) -DARB_USE_BUNDLED_LIBS=ON -DARB_WITH_PYTHON=ON ..
# make -j 12 && make install

import sys
import importlib.util
import os.path

#### IF INSTALLED USING PIP:
ARBOR_LOCATION = 'default'
ARBOR_BUILD_CATALOGUE = 'arbor-build-catalogue'

#### IF INSTALLED FROM SOURCE:
#ARBOR_REPOSITORY = '/home/llandsmeer/tmp/a-e/arbor'
#ARBOR_LOCATION = os.path.join(ARBOR_REPOSITORY,
#        'build/prefix/lib/python3.10/site-packages/arbor/__init__.py')
#ARBOR_BUILD_CATALOGUE = os.path.join(ARBOR_REPOSITORY,
#        'build/prefix/bin/arbor-build-catalogue')

if ARBOR_LOCATION != 'default':
    spec = importlib.util.spec_from_file_location('arbor', ARBOR_LOCATION)
    module = importlib.util.module_from_spec(spec)
    sys.modules[spec.name] = module
    spec.loader.exec_module(module)
# END INSTALLATION OPTIONS

import numpy as np
import re
import subprocess
import tempfile
import arbor
import matplotlib.pyplot as plt

mechanism_source = '''
NEURON {
  SUFFIX noise
  NONSPECIFIC_CURRENT i
  RANGE mu, theta, sigma
}

PARAMETER {
  mu = 0
  theta = 1
  sigma = 1
}

STATE { x }

WHITE_NOISE { W }

DERIVATIVE dX {
}

BREAKPOINT {
  SOLVE dX METHOD stochastic
  i = W
}
'''

def build_noise_catalogue():
    with tempfile.TemporaryDirectory() as tmpdir:
        with open(os.path.join(tmpdir, 'noise.mod'), 'w') as f:
            print(mechanism_source, file=f)
        out = subprocess.getoutput(f'{ARBOR_BUILD_CATALOGUE} noise {tmpdir}')
        m = re.search('[^ ]*\.so', out)
        if m is None:
            print(out)
            input()
            exit(1)
        return arbor.load_catalogue(m.group(0))


def make_cable_cell(gid):
    tree = arbor.segment_tree()
    soma = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
    labels = arbor.label_dict(dict(soma="(tag 1)", root= "(root)"))
    decor = arbor.decor()
    decor.paint('"soma"', arbor.density("noise"))
    return arbor.cable_cell(tree, labels, decor)

class Recipe(arbor.recipe):
    def __init__(self, ncells=100):
        arbor.recipe.__init__(self)
        self.ncells = ncells
        self.props = arbor.neuron_cable_properties()
        self.props.catalogue.extend(build_noise_catalogue(), '')
    def num_cells(self): return self.ncells
    def cell_description(self, gid): return make_cable_cell(gid)
    def cell_kind(self, gid): return arbor.cell_kind.cable
    def connections_on(self, gid): return []
    def event_generators(self, gid): return []
    def probes(self, gid): return [arbor.cable_probe_membrane_voltage('"root"')]
    def global_properties(self, kind): return self.props

# (11) Instantiate recipe
recipe = Recipe(ncells=100)
ctx = arbor.context('avail_threads')
sim = arbor.simulation(recipe, ctx)
handles = [sim.sample((gid, 0), arbor.regular_schedule(1)) for gid in range(recipe.num_cells())]
sim.run(100)
voltages = np.array([sim.samples(handles[gid])[0][0][:, 1] for gid in range(recipe.num_cells())])

plt.plot(voltages.T)
plt.show()