Changeset 955

Timestamp:

04/22/11 18:04:21 (13 months ago)

Author:

apdavison

Message:

Cell types specified in 9ML can now have multiple synaptic channels (pyNN.neuron)

Location:

trunk

Files:

• examples/nineml_neuron.py (modified) (6 diffs)
• src/connectors.py (modified) (1 diff)
• src/neuron/nineml.py (modified) (12 diffs)
• src/nineml/read.py (modified) (2 diffs)

trunk/examples/nineml_neuron.py

@@ +1 @@
+"""
+Example of using a cell type defined in 9ML with pyNN.neuron
+"""
+
+
 import sys
-sys.path.append("/home/andrew/dev/nineml_all/trunk/lib9ml/python/examples/AL/")
-sys.path.append("/home/andrew/dev/nineml_all/trunk/code_generation/nmodl/")
+from os.path import abspath, realpath, join
+import nineml
+root = abspath(join(realpath(nineml.__path__[0]), "../../.."))
+sys.path.append(join(root, "lib9ml/python/examples/AL"))
+sys.path.append(join(root, "code_generation/nmodl"))                
 leaky_iaf = __import__("leaky_iaf")
 coba_synapse = __import__("coba_synapse")
@@ -7 +15 @@
 from pyNN.utility import init_logging
 import pyNN.neuron as sim
+from copy import deepcopy
 
 init_logging(None, debug=True)
@@ -14 +23 @@
                                 leaky_iaf.c1,
                                 excitatory=coba_synapse.c1,
-                                port_map=[('V', 'V'), ('Isyn', 'Isyn')])
+                                inhibitory=deepcopy(coba_synapse.c1),
+                                port_map={
+                                    'excitatory': [('V', 'V'), ('Isyn', 'Isyn')],
+                                    'inhibitory': [('V', 'V'), ('Isyn', 'Isyn')]
+                                },
+                                weight_variables={
+                                    'excitatory': 'q',
+                                    'inhibitory': 'q'
+                                })
 
 parameters = {
@@ -20 +37 @@
     'gL': 50.0,
     't_ref': 5.0,
-    'tau': 2.0,
+    'excitatory_tau': 2.0,
+    'inhibitory_tau': 5.0,
+    'excitatory_E': 0.0,
+    'inhibitory_E': -70.0,
     'theta': -50.0,
     'vL': -65.0,
@@ -26 +46 @@
 }
 
-#celltype = celltype_cls(parameters)
-#cell = celltype.model(**celltype.parameters)
-
 cells = sim.Population(1, celltype_cls, parameters)
 cells.initialize('V', parameters['vL'])
 cells.initialize('t_spike', -1e99) # neuron not refractory at start
+cells.initialize('regime', 1002) # temporary hack
 
-input = sim.Population(1, sim.SpikeSourcePoisson, {'rate': 100})
+input = sim.Population(2, sim.SpikeSourcePoisson, {'rate': 100})
 
-connector = sim.OneToOneConnector(weights=0.1, delays=0.5)
-conn = sim.Projection(input, cells, connector, target='excitatory')
+connector = sim.OneToOneConnector(weights=1.0, delays=0.5)
+conn = [sim.Projection(input[0:1], cells, connector, target='excitatory'),
+        sim.Projection(input[1:2], cells, connector, target='inhibitory')]
 
 cells._record('V')
-cells._record('g')
+cells._record('excitatory_g')
+cells._record('inhibitory_g')
 cells.record()
 
@@ -45 +65 @@
 
 cells.recorders['V'].write("Results/nineml_neuron.V", filter=[cells[0]])
-cells.recorders['g'].write("Results/nineml_neuron.g", filter=[cells[0]])
+cells.recorders['excitatory_g'].write("Results/nineml_neuron.g_exc", filter=[cells[0]])
+cells.recorders['inhibitory_g'].write("Results/nineml_neuron.g_inh", filter=[cells[0]])
 
 sim.end()

trunk/src/connectors.py

@@ -442 +442 @@
         `weights` -- may either be a float, a RandomDistribution object, a list/
                      1D array with at least as many items as connections to be
-                     created, or a DistanceDependence object. Units nA.
+                     created, or a distance expression as for `d_expression`. Units nA.
         `delays`  -- as `weights`. If `None`, all synaptic delays will be set
                      to the global minimum delay.

trunk/src/neuron/nineml.py

@@ -1 @@
-## Not compatible with Python 2.4 ###
-from __future__ import absolute_import
+"""
+Support cell types defined in 9ML with NEURON.
+
+Requires the 9ml2nmodl script to be on the path.
+
+Classes:
+    NineMLCell       - a single neuron instance
+    NineMLCellType   - base class for cell types, not used directly
+
+Functions:
+    nineml_cell_type - return a new NineMLCellType subclass
+
+Constants:
+    NMODL_DIR        - subdirectory to which NMODL mechanisms will be written
+    
+"""
+
+from __future__ import absolute_import # Not compatible with Python 2.4
 import subprocess
 import neuron
@@ -7 +23 @@
 import logging
 import os
+import re
 from itertools import chain
 
@@ -60 +77 @@
 
 
-def _compile_nmodl(nineml_component):
+def _compile_nmodl(nineml_component, weight_variables): # weight variables should really be within component
+    """
+    Generate NMODL code for the 9ML component, run "nrnivmodl" and then load
+    the mechanisms into NEURON.
+    """
     if not os.path.exists(NMODL_DIR):
         os.makedirs(NMODL_DIR)
@@ -69 +90 @@
     nineml_component.write(xml_file)
     nineml2nmodl = __import__("9ml2nmodl")
-    nineml2nmodl.write_nmodl(xml_file)
+    nineml2nmodl.write_nmodl(xml_file, weight_variables) # weight variables should really come from xml file
     p = subprocess.check_call(["nrnivmodl"])
     os.chdir(cwd)
@@ -75 +96 @@
 
 
+def _add_prefix(synapse_model, prefix, port_map):
+    """
+    Add a prefix to all variables in `synapse_model`, except for variables with
+    receive ports and specified in `port_map`.
+    """
+    synapse_model.__cache__ = {}
+    exclude = []
+    new_port_map = []
+    for name1, name2 in port_map:
+        if synapse_model.ports_map[name2].mode == 'recv':
+            exclude.append(name2)
+            new_port_map.append((name1, name2))
+        else:
+            new_port_map.append((name1, prefix + '_' + name2))
+    synapse_model.add_prefix(prefix + '_', exclude=exclude)
+    return new_port_map
+
+
 class _build_nineml_celltype(type):
     """
@@ -80 +119 @@
     """
     def __new__(cls, name, bases, dct):
-        assert len(dct["synapse_models"]) == 1, "For now, can't handle multiple synapse models"
-        combined_model = join(dct["neuron_model"],
-                              dct["synapse_models"].values()[0],
-                              dct["port_map"],
-                              name=name)
+        # join the neuron and synapse components into a single component
+        combined_model = dct["neuron_model"]
+        for label in dct["synapse_models"].keys():
+            port_map = dct["port_map"][label]
+            port_map = _add_prefix(dct["synapse_models"][label], label, port_map)
+            dct["weight_variables"][label] = label + "_" + dct["weight_variables"][label]
+            combined_model = join(combined_model,
+                                  dct["synapse_models"][label],
+                                  port_map,
+                                  name=name)
         dct["combined_model"] = combined_model
+        # set class attributes required for a PyNN cell type class
         dct["default_parameters"] = dict((name, 1.0)
                                       for name in combined_model.parameters)
@@ -92 +137 @@
         dct["synapse_types"] = dct["synapse_models"].keys() #really need an ordered dict
         dct["injectable"] = True # need to determine this. How??
-        dct["recordable"] = [port.name for port in combined_model.analog_ports] + ['spikes']
+        dct["recordable"] = [port.name for port in combined_model.analog_ports] + ['spikes', 'regime']
         dct["standard_receptor_type"] = (dct["synapse_types"] == ('excitatory', 'inhibitory'))
         dct["conductance_based"] = True # how to determine this??
         dct["model_name"] = name
         logger.debug("Creating class '%s' with bases %s and dictionary %s" % (name, bases, dct))
-        _compile_nmodl(combined_model)
+        # generate and compile NMODL code, then load the mechanism into NEUORN
+        _compile_nmodl(combined_model, dct["weight_variables"]) # weight variables should really be stored within combined_model
         return type.__new__(cls, name, bases, dct)
     
     
 
-def nineml_cell_type(name, neuron_model, port_map={}, **synapse_models):
+def nineml_cell_type(name, neuron_model, port_map={}, weight_variables={}, **synapse_models):
     """
     Return a new NineMLCellType subclass.
@@ -109 +155 @@
                                   {'neuron_model': neuron_model,
                                    'synapse_models': synapse_models,
-                                   'port_map': port_map})
-    
+                                   'port_map': port_map,
+                                   'weight_variables': weight_variables})
+
 
 def join(c1, c2, port_map=[], name=None):
@@ -116 +163 @@
     logger.debug("Joining components %s and %s with port map %s" % (c1, c2, port_map))
     logger.debug("New component will have name '%s'" % name)
-    bindings = [] # TODO: combine bindings from c1 and c2
+    # combine bindings from c1 and c2
+    bindings = {}
+    for b in chain(c1.bindings, c2.bindings):
+        bindings[b.name] = b
+    # combine ports (some will later be removed)
     all_ports = c1.ports_map.copy()
     all_ports.update(c2.ports_map)
+    # event ports do not be passed to the constructor, as they are attached to transitions
     for port_name, port in all_ports.items():
         if isinstance(port, nineml.EventPort):
             all_ports.pop(port_name)
+    # connect ports.
+    # currently, when ports are connected they disappear. It might be better to
+    # explicitly keep the ports in the new component but mark them as connected
     for name1, name2 in port_map:
-        assert name1 in c1.ports_map
-        assert name2 in c2.ports_map
-        all_ports.pop(name1)
-        if name1 != name2:
-            #c2.substitute(name2, name1) # need to implement this
-            all_ports.pop(name2)
-        
+        assert name1 in c1.ports_map, "%s is not in %s" % (name1, c1.ports_map.keys())
+        assert name2 in c2.ports_map, "%s is not in %s" % (name2, c2.ports_map.keys())
+
         port1 = c1.ports_map[name1]
         port2 = c2.ports_map[name2]
@@ -135 +186 @@
         if port1.mode == 'send':
             send_port = port1
-            port_name = name1
+            recv_port = port2
+            send_port_name = name1
+            recv_port_name = name2
         else:
             send_port = port2
-            port_name = name2
+            recv_port = port1
+            send_port_name = name2
+            recv_port_name = name1
+        # when connecting ports in which the send port has an expression, need
+        # to create a binding for this expression in the new component
         if send_port.expr:
             func_args = c1.non_parameter_symbols.union(c2.non_parameter_symbols).intersection(send_port.expr.names)
-            lhs = "%s(%s)" % (port_name, ",".join(func_args))
-            bindings.append(nineml.Binding(lhs, send_port.expr.rhs))
+            lhs = "%s(%s)" % (send_port_name, ",".join(func_args))
+            send_binding = nineml.Binding(lhs, send_port.expr.rhs)
+            bindings[send_binding.name] = send_binding
             for eq in chain(c1.equations, c2.equations):
-                if port_name in eq.names:
-                    eq.rhs = eq.rhs_name_transform({port_name: lhs})
+                if send_port_name in eq.names:
+                    eq.rhs = eq.rhs_name_transform({send_port_name: lhs})            
+            if recv_port.mode == 'reduce':
+                # need to retain reduce ports as they can be connected to in a future join
+                if recv_port_name in bindings:
+                    # this reduce port has already been connected to, so combine using its reduce_op
+                    reduce_binding = bindings[recv_port_name]
+                    func_args = func_args.union(reduce_binding.args)
+                    lhs = "%s(%s)" % (recv_port_name, ",".join(func_args))
+                    rhs = recv_port.reduce_op.join([reduce_binding.rhs, send_binding.lhs])
+                else:
+                    # this is the first time this reduce port has been connected to
+                    lhs = "%s(%s)" % (recv_port_name, ",".join(func_args))
+                    rhs = send_binding.lhs
+                bindings[recv_port_name] = nineml.Binding(lhs, rhs)
+                recv_port.connected = True
+            else:
+                all_ports.pop(name1)
+        else:
+            if recv_port.mode == 'reduce':
+                raise NotImplementedError
+            else:
+                all_ports.pop(name1)
+
+        if name1 != name2:
+            #c2.substitute(name2, name1) # need to implement this. Currently this all only works if name1 == name2
+                                         # probably needs to happen sooner in the function
+            all_ports.pop(name2)
+
+    # where parameters have become bindings due to connecting ports, replace
+    # bare names with function calls in the equations
+    for bname, binding in bindings.items():
+        for eq in chain(c1.equations, c2.equations):
+            if bname in eq.names:
+                print "#### replacing %s by %s" % (bname, binding.lhs)
+                pattern = re.compile(r'%s(\([\w\, ]*\))?' % bname)
+                m = pattern.search(eq.rhs)
+                if m:
+                    eq.rhs = pattern.sub(binding.lhs, eq.rhs)
+                else:
+                    eq.rhs = eq.rhs_name_transform({bname: binding.lhs})
+
+    # create new regimes from all possible combinations of the regimes from the
+    # two components
     regime_map = {}
     for r1 in c1.regimes:
         regime_map[r1.name] = {}
         for r2 in c2.regimes:
-            kwargs = {'name': "%s_AND_%s" % (r1.name, r2.name)}
+            if r1.name == r2.name:
+                new_name = r1.name
+            else:
+                new_name = "%s_AND_%s" % (r1.name, r2.name)
+            kwargs = {'name': new_name}
             new_regime = nineml.Regime(*r1.nodes.union(r2.nodes), **kwargs)
             regime_map[r1.name][r2.name] = new_regime
+    # create transitions between all the new regimes
     transitions = []
     for r1 in c1.regimes:
@@ -168 +273 @@
                                                    condition=t.condition)
                 transitions.append(new_transition)
+
     regimes = []
     for d in regime_map.values():
@@ -176 +282 @@
                             transitions=transitions,
                             ports=all_ports.values(),
-                            bindings=bindings)
+                            bindings=bindings.values())

trunk/src/nineml/read.py

@@ -188 +188 @@
     """
     
-    def __init__(self, sim, nineml_file):
+    def __init__(self, sim, nineml_model):
         """
         Instantiate a network from a 9ML file, in the specified simulator.
         """
         global random_distributions
-        assert isinstance(nineml_file, basestring)
         self.sim = sim
-        self.nineml_model = nineml.parse(nineml_file)
+        if isinstance(nineml_model, basestring):
+            self.nineml_model = nineml.parse(nineml_model)
+        elif isinstance(nineml_model, nineml.Model):
+            self.nineml_model = nineml_model
+        else:
+            raise TypeError("nineml_model must be a nineml.Model instance or the path to a NineML XML file.")
         self.random_distributions = {}
         self.assemblies = {}
@@ -207 +211 @@
     def _handle_group(self, group):
         # create an Assembly
-        self.assemblies[group.name] = self.sim.Assembly(group.name)
+        self.assemblies[group.name] = self.sim.Assembly(label=group.name)
         
         # extract post-synaptic response definitions from projections