Update/workshop 2024

bmtk/builder/builder_utils.py

@@ -87,6 +87,8 @@ def check_properties_across_ranks(properties, graph_type='node'):
                 phash = hashlib.md5(pval.encode('utf-8')).hexdigest()
 
             elif isinstance(pval, (int, float, bool)):
+                if np.isnan(pval):
+                    pval = 'NONE'
                 phash = pval
 
             elif isinstance(pval, (list, tuple)):

bmtk/simulator/bionet/biosimulator.py

@@ -376,41 +376,8 @@ def from_config(cls, config, network, set_recordings=True):
                 sim.add_mod(mods.IClampMod(input_type=sim_input.input_type, **sim_input.params))
 
             elif sim_input.module == "SEClamp":
-                node_set = network.get_node_set(sim_input.node_set)
-                try:
-                    len(sim_input.params['amps'])
-                except:
-                    sim_input.params['amps']=[float(sim_input.params['amps'])]
+                sim.add_mod(mods.SEClamp(input_type=sim_input.input_type, **sim_input.params))
 
-                try: 
-                    len(sim_input.params['durations'])
-                except:
-                    sim_input.params['durations']=[float(sim_input.params['durations'])]
-
-                amplitudes = sim_input.params['amps']
-                durations = sim_input.params['durations']
-                rs = None
-
-                if "rs" in sim_input.params.keys():
-                    try: 
-                        len(sim_input.params['rs'])
-                    except:
-                        sim_input.params['rs']=[float(sim_input.params['rs'])]
-                    if len(sim_input.params['rs'])>1:
-                        sim_input.params['rs']=[float(i) for i in sim_input.params['rs']]
-                    rs = sim_input.params["rs"]
-
-                try:
-                    sim_input.params['gids']
-                except:
-                    sim_input.params['gids'] = None
-                if sim_input.params['gids'] is not None:
-                    gids = sim_input.params['gids']
-                else:
-                    gids = list(node_set.gids())
-
-                sim.attach_se_voltage_clamp(amplitudes, durations, gids, rs)
-
             elif sim_input.module == 'xstim':
                 sim.add_mod(mods.XStimMod(**sim_input.params))
 

bmtk/simulator/bionet/modules/__init__.py

@@ -30,3 +30,4 @@
 from .iclamp import IClampMod
 from .comsol import ComsolMod
 from .ecephys_module import BioECEphysUnitsModule
+from .seclamp import SEClamp

bmtk/simulator/bionet/modules/seclamp.py

@@ -0,0 +1,70 @@
+import numpy as np
+from neuron import h
+
+from bmtk.simulator.bionet.modules.sim_module import SimulatorMod
+from bmtk.simulator.bionet.io_tools import io
+from bmtk.simulator.core.modules import iclamp
+
+
+class SEClamp(SimulatorMod):
+    def __init__(self, input_type, **mod_args):
+
+        # Select location to place iclamp, if not specified use the center of the soma
+        self._section_name = mod_args.get('section_name', 'soma')
+        self._section_index = mod_args.get('section_index', 0)
+        self._section_dist = mod_args.get('section_dist', 0.5)
+
+        self._rs = mod_args.get('resistance', 1.0)
+        self._vc = mod_args.get('vc', 0.0)
+
+        self._delay = mod_args['delay']
+        self._amp = mod_args['amp']
+        self._duration = mod_args['duration']
+        self._ton = self._delay
+        self._toff = self._delay + self._duration
+
+        # Check f section_index is a range (ie. "section_index": [500.0, 1000.0])
+        self._ranged_index = isinstance(self._section_index, (list, tuple))
+
+        # SEClamp objects need to be saved in memory otherwise NRN will try to garbage collect them
+        # prematurly
+        self._seclamp = None
+
+        self._node_set = mod_args.get('node_set', 'all') 
+
+    def initialize(self, sim):
+        select_gids = list(sim.net.get_node_set(self._node_set).gids())
+        gids_on_rank = list(set(select_gids) & set(select_gids))
+
+        for gid in gids_on_rank:           
+            cell = sim.net.get_cell_gid(gid)
+
+            if self._ranged_index:
+                hobj_sec = self._find_section(cell)
+            else:
+                hobj_sec = getattr(cell.hobj, self._section_name)[self._section_index](self._section_dist)
+
+            self._seclamp = self.create_clamp(hobj_sec)
+
+    def step(self, sim, tstep):
+        if self._ton <= sim.dt*tstep <= self._toff:
+            self._seclamp.rs = self._rs
+            # print(sim.dt*tstep)
+        else:
+            self._seclamp.rs = 10.0e20
+
+
+    def create_clamp(self, hobj):
+        stim = h.SEClamp(hobj)
+
+        stim.dur1 = self._delay
+        stim.dur2 = self._duration
+        stim.dur3 = 0.0
+
+        stim.amp1 = 0.0
+        stim.amp2 = self._amp
+        stim.amp3 = 0.0
+        stim.vc = self._vc
+
+        return stim
+

bmtk/simulator/bionet/modules/xstim_waveforms.py

@@ -56,7 +56,7 @@ def calculate(self, t):  # TODO better name
 class WaveformCustom(BaseWaveform):
     """Custom waveform defined by csv file"""
     def __init__(self, waveform_file):
-        self.definition = pd.read_csv(waveform_file, sep='\t')
+        self.definition = pd.read_csv(waveform_file, sep=' ')
 
     def calculate(self, t):
         return np.interp(t, self.definition["time"], self.definition["amplitude"])

bmtk/simulator/filternet/filtersimulator.py

@@ -48,6 +48,16 @@ def add_movie(self, movie_type, params):
                 else:
                     raise Exception('Could not find movie "data_file" in config to use as input.')
 
+                # If file passed in is a npz compressed file then it is in a dictionary format and need to find
+                # the key-value pair containing the array
+                if isinstance(m_data, np.lib.npyio.NpzFile):
+                    try:
+                        for key in m_data:
+                            m_data = m_data[key]
+                            break
+                    except IndexError as ie:
+                        io.log_warning('Was unable to find array from compressed numpy matrix file.')
+
                 # contrast_min, contrast_max = m_data.min(), m_data.max()
                 normalize_data = params.get('normalize', False)
                 if normalize_data:
@@ -181,21 +191,21 @@ def run(self):
         ten_percent = int(np.ceil(n_cells_on_rank*0.1))
         rank_msg = '' if bmtk_world_comm.MPI_size < 2 else ' (on rank {})'.format(bmtk_world_comm.MPI_rank)
 
-        max_fr = np.empty(len(cells_on_rank))
+        # max_fr = np.empty(len(cells_on_rank))
         for cell_num, cell in enumerate(cells_on_rank):
             for movie, options in zip(self._movies, self._eval_options):
                 if cell_num > 0 and cell_num % ten_percent == 0:
                     io.log_debug(' Processing cell {} of {}{}.'.format(cell_num, n_cells_on_rank, rank_msg))
                 ts, f_rates = cell.lgn_cell_obj.evaluate(movie, **options)
-                max_fr[cell_num] = np.max(f_rates)
+                # max_fr[cell_num] = np.max(f_rates)
                 if movie.padding:
                     f_rates = f_rates[int(movie.data.shape[0]-movie.data_orig.shape[0]):]
                     ts = ts[int(movie.data.shape[0]-movie.data_orig.shape[0]):]
                     ts = ts-ts[0]
 
                 for mod in self._sim_mods:
                     mod.save(self, cell, ts, f_rates)
-        io.log_info('Max firing rate: {}'.format(np.max(max_fr)))
+        # io.log_info('Max firing rate: {}'.format(np.max(max_fr)))
         io.log_info('Done.')
         for mod in self._sim_mods:
             mod.finalize(self)

bmtk/simulator/filternet/lgnmodel/cursor.py

@@ -225,14 +225,20 @@ def evaluate(self, threshold=0):
         full_temporal_kernel = self.temporal_kernel.full()
 
         # Convolve every frame in the movie with the spatial filter. First find the range of rows (range min and max)
-        #  and columns in the filter that are above threshold, that way only portion of movie/filter are multiplied
+        # and columns in the filter that are above threshold, that way only portion of movie/filter are multiplied
         # together
 
         # Using > instead of >= makes the code faster if there are many zeros in the
         # spatial kernel. This will not affect the results.
         nonzero_inds = np.where(np.abs(full_spatial_kernel[0, :, :]) > threshold)
-        rm, rM = nonzero_inds[0].min(), nonzero_inds[0].max()
-        cm, cM = nonzero_inds[1].min(), nonzero_inds[1].max()
+        if len(nonzero_inds[0]) == 0:
+            # If spatial kernel is all 0's then don't try to segment the receptive filed. Use full frame during convolve.
+            rm, rM = 0, 0
+            cm, cM = full_spatial_kernel.shape[1], full_spatial_kernel.shape[2]
+        else:
+            rm, rM = nonzero_inds[0].min(), nonzero_inds[0].max()
+            cm, cM = nonzero_inds[1].min(), nonzero_inds[1].max()
+
         convolution_answer_sep_spatial = (self.movie.data[:, rm:rM+1, cm:cM+1] *
                                           full_spatial_kernel[:, rm:rM+1, cm:cM+1]).sum(axis=1).sum(axis=1)
 

bmtk/simulator/pointnet/modules/spikes_inputs.py

@@ -24,7 +24,7 @@
 from bmtk.simulator.pointnet.modules.sim_module import SimulatorMod
 from bmtk.simulator.pointnet.io_tools import io
 from bmtk.utils.reports.spike_trains import SpikeTrains
-from bmtk.simulator.bionet.pyfunction_cache import py_modules
+from bmtk.simulator.pointnet.pyfunction_cache import py_modules
 
 
 class SpikesInputsMod(SimulatorMod):
@@ -48,7 +48,7 @@ def initialize(self, sim):
                 io.log_exception(f'Could not find @spikes_generator function "{spikes_generator}" required for {self._name} inputs.')
             spikes_func = py_modules.spikes_generator(name=spikes_generator)
 
-            self._spike_trains = SpikeTrains(cace_to_disk=False)
+            self._spike_trains = SpikeTrains(cache_to_disk=False)
             for node in node_set.fetch_nodes():
                 timestamps = spikes_func(node, sim)
                 self._spike_trains.add_spikes(

bmtk/simulator/pointnet/pointnetwork.py

@@ -86,6 +86,8 @@ def __init__(self, **properties):
         self._nodes_table = {}
         self._gid2nestid = {}
 
+        self._nest_modules = []
+
         self._gid_map = GidPool()
         self._virtual_gids = GidPool()
 
@@ -131,6 +133,10 @@ def add_weight_function(self, fnc, name=None, **kwargs):
         fnc_name = name if name is not None else function.__name__
         self.__weight_functions[fnc_name] = functools.partial(fnc)
 
+    def add_nest_module(self, module_path):
+        if module_path not in self._nest_modules:
+            self._nest_modules.append(module_path)
+
     def set_default_weight_function(self, fnc):
         self.add_weight_function(fnc, 'default_weight_fnc', overwrite=True)
 

bmtk/simulator/pointnet/pointsimulator.py

@@ -25,6 +25,8 @@
 import nest
 from six import string_types
 from six import moves
+from pathlib import Path
+import platform
 
 from bmtk.simulator.core.simulator import Simulator
 from bmtk.simulator.pointnet.config import Config
@@ -72,6 +74,7 @@ def __init__(self, graph, dt=0.001, overwrite=True, print_time=False, n_thread=1
         # TODO: move this into it's own function and make sure it is called before network is built
         nest.ResetKernel()
         nest.SetKernelStatus({"resolution": self._dt, "overwrite_files": self._overwrite, "print_time": print_time, "local_num_threads": n_thread})
+        self._load_nest_modules()
 
     @property
     def tstart(self):
@@ -131,6 +134,54 @@ def _get_block_trial(self, duration):
             data_res = -1
         return n, res, data_res
 
+    def _add_library_path(self, lib_path):
+        if isinstance(lib_path, Path):
+            lib_path = lib_path.as_posix()
+        system = platform.system()
+        env_var = "LD_LIBRARY_PATH" if system in ['Linux', 'Windows', 'Java'] else 'DYLD_LIBRARY_PATH'
+        env_val = os.environ.get(env_var, '')
+
+        if lib_path not in env_val:
+            os.environ[env_var] = os.pathsep.join([lib_path, env_val])
+
+
+    def _load_nest_modules(self):
+        # If there is a "nest_modules" entry in the configuration "components" then go through and add them to
+        # the network _nest_modules list. They will be processed the same as calling network.add_nest_module().
+        components_path = self.net._components.get('nest_modules', [])
+        if isinstance(components_path, str):
+            self.net.add_nest_module(components_path)
+        elif isinstance(components_path, list):
+            for cpath in components_path:
+                self.net.add_nest_module(cpath)
+        else:
+            raise ValueError('Unable to load components/nest_modules value.')
+
+        # Go through all added nest modules and try to call nest.Install() for them. 
+        for module in self.net._nest_modules:
+            if Path(module).is_dir():
+                # If module is a directory then add path to LD_LIBRARY_PATH then try to load all the .so/.a library binaries in directory
+                lib_dir = Path(module).resolve()
+                self._add_library_path(lib_dir)
+                so_search = Path(lib_dir) / '*.so'
+                a_search = Path(lib_dir) / '*.a'
+
+                for lib_file in glob.glob(so_search.as_posix()) + glob.glob(a_search.as_posix()): # lib_path.parent.resolve()
+                    module_name = Path(lib_file).name
+                    nest.Install(module_name)
+
+            elif Path(module).is_file():
+                # If user tries to pass in a path to a library binary
+                lib_path = Path(module)
+                lib_dir = lib_path.parent.resolve()
+                lib_filename = lib_path.name
+                self._add_library_path(lib_dir)
+                nest.Install(lib_filename)               
+
+            else:
+                # If user just tries nest.Install('mymodule')
+                nest.Install(module)
+
     '''
     def set_spikes_recordings(self):
         # TODO: Pass in output-dir and file name to save to

bmtk/utils/create_environment/__main__.py

@@ -52,7 +52,7 @@ def __split_list(options, opt_name):
     parser.add_option('-c', '--components-dir', dest='components_dir', default=None,
                       help="Directory to use for parameter files, morphology, and other used components.")
     parser.add_option('--tstop', type='float', dest='tstop', default=1000.0)
-    parser.add_option('--dt', type=float, dest='dt', help='simulation time step dt', default=0.001)
+    parser.add_option('--dt', type=float, dest='dt', help='simulation time step dt', default=None)
     parser.add_option('--report-vars', dest='report_vars', type='string', action='callback',
                       callback=__list_parser, default=[],
                       help='A list of membrane variables to record from; v, cai, etc.')

bmtk/utils/create_environment/create_environment.py

@@ -41,7 +41,7 @@ def create_environment(simulator,
                        se_voltage_clamp=None,
                        tstart=0.0,
                        tstop=1000.0,
-                       dt=0.001,
+                       dt=None,
                        dL=20.0,
                        spikes_threshold=-15.0,
                        nsteps_block=5000,
@@ -50,7 +50,7 @@ def create_environment(simulator,
                        compile_mechanisms=False,
                        use_relative_paths=True,
                        include_examples=False
-                       ):
+                      ):
     """
 
     :param simulator:

bmtk/utils/create_environment/env_builder.py

@@ -776,6 +776,15 @@ def target_simulator(self):
     @property
     def bmtk_simulator(self):
         return 'filternet'
+
+    def _add_run_params(self, tstart=0.0, tstop=1000.0, dt=None, **kwargs):
+        self._simulation_config['run'] = {
+            'tstart': tstart,
+            'tstop': tstop
+        }
+
+        if dt is not None:
+            self._simulation_config['run']['dt'] = dt
 
     def _add_output_section(self):
         super(FilterNetEnvBuilder, self)._add_output_section()

examples/point_nestml_izh/config.circuit.json

@@ -0,0 +1,38 @@
+{
+  "manifest": {
+	"$BASE_DIR": ".",
+    "$NETWORK_DIR": "$BASE_DIR/network",
+    "$MODELS_DIR": "$BASE_DIR/../point_components"
+  },
+
+  "components": {
+    "point_neuron_models_dir": "$MODELS_DIR/cell_models",
+	  "synaptic_models_dir": "$MODELS_DIR/synaptic_models",
+    "nest_modules": "components/nestml/nestml_izh_module.so"
+  },
+
+  "networks": {
+    "nodes": [
+      {
+        "nodes_file": "$NETWORK_DIR/cortex_nodes.h5",
+        "node_types_file": "$NETWORK_DIR/cortex_node_types.csv"
+      },
+      {
+        "nodes_file": "$NETWORK_DIR/thalamus_nodes.h5",
+        "node_types_file": "$NETWORK_DIR/thalamus_node_types.csv"
+      }
+    ],
+    "edges": [
+      {
+        "edges_file": "$NETWORK_DIR/cortex_cortex_edges.h5",
+        "edge_types_file": "$NETWORK_DIR/cortex_cortex_edge_types.csv",
+        "enabled": false
+      },
+      {
+        "edges_file": "$NETWORK_DIR/thalamus_cortex_edges.h5",
+        "edge_types_file": "$NETWORK_DIR/thalamus_cortex_edge_types.csv",
+        "enabled": true
+      }
+    ]
+  }
+}