from .rxdmath import _Arithmeticed import weakref from .section1d import Section1D from neuron import h, nrn from . import node, nodelist, rxdmath, region import numpy import warnings import itertools from .rxdException import RxDException from . import initializer import collections # The difference here is that defined species only exists after rxd initialization _all_species = [] _defined_species = {} def _get_all_species(): return _defined_species _species_count = 0 _has_1d = False _has_3d = False def _1d_submatrix_n(): if not _has_1d: return 0 elif not _has_3d: return len(node._states) else: return numpy.min([sp()._indices3d() for sp in list(_get_all_species().values()) if sp() is not None]) class _SpeciesMathable(object): # support arithmeticing def __neg__(self): return -1 * _Arithmeticed(self) def __Pos__(self): return _Arithmeticed(self) def __pow__(self, other): return _Arithmeticed(self) ** other def __add__(self, other): return _Arithmeticed(self) + other def __sub__(self, other): return _Arithmeticed(self) - other def __mul__(self, other): return _Arithmeticed(self) * other def __div__(self, other): return _Arithmeticed(self) / other def __radd__(self, other): return _Arithmeticed(other) + self def __rmul__(self, other): return _Arithmeticed(self) * other def __rdiv__(self, other): return _Arithmeticed(other) / self def __rsub__(self, other): return _Arithmeticed(other) - self def __abs__(self): return abs(_Arithmeticed(self)) def _evaluate(self, location): return _Arithmeticed(self)._evaluate(location) def __ne__(self, other): other = rxdmath._ensure_arithmeticed(other) self2 = rxdmath._ensure_arithmeticed(self) rxdmath._validate_reaction_terms(self2, other) return rxdmath._Reaction(self2, other, '<>') def __gt__(self, other): other = rxdmath._ensure_arithmeticed(other) self2 = rxdmath._ensure_arithmeticed(self) rxdmath._validate_reaction_terms(self2, other) return rxdmath._Reaction(self2, other, '>') def __lt__(self, other): other = rxdmath._ensure_arithmeticed(other) self2 = rxdmath._ensure_arithmeticed(self) rxdmath._validate_reaction_terms(self2, other) return rxdmath._Reaction(self2, other, '<') @property def _semi_compile(self): return 'species%d' % (self._id) def _involved_species(self, the_dict): the_dict[self._semi_compile] = weakref.ref(self) @property def d(self): """diffusion constant. write-only""" raise RxDException('diffusion constant is write-only') @d.setter def d(self, value): from . import rxd if hasattr(self, '_allow_setting'): self._d = value else: _volumes, _surface_area, _diffs = node._get_data() _diffs[self.indices()] = value rxd._setup_matrices() class SpeciesOnRegion(_SpeciesMathable): def __init__(self, species, region): """The restriction of a Species to a Region.""" global _species_count self._species = weakref.ref(species) self._region = weakref.ref(region) self._id = _species_count _species_count += 1 def __eq__(self, other): """test for equality. Two SpeciesOnRegion objects are equal if they refer to the same species on the same region and both the species and the region still exist. """ return (self._species() == other._species()) and (self._region() == other._region()) and (self._species() is not None) and (self._region() is not None) def __hash__(self): # TODO: replace this to reduce collision risk; how much of a problme is that? return 1000 * (hash(self._species()) % 1000) + (hash(self._region()) % 1000) def __repr__(self): return '%r[%r]' % (self._species(), self._region()) def _short_repr(self): return '%s[%s]' % (self._species()._short_repr(), self._region()._short_repr()) def __str__(self): return '%s[%s]' % (self._species()._short_repr(), self._region()._short_repr()) def indices(self, r=None, secs=None): """If no Region is specified or if r is the Region specified in the constructor, returns a list of the indices of state variables corresponding to the Species when restricted to the Region defined in the constructor. If r is a different Region, then returns the empty list. If secs is a set return all indices on the regions for those sections. """ #if r is not None and r != self._region(): # raise RxDException('attempt to access indices on the wrong region') # TODO: add a mechanism to catch if not right region (but beware of reactions crossing regions) if self._species() is None or self._region() is None: return [] else: return self._species().indices(self._region(), secs) def __getitem__(self, r): """Return a reference to those members of this species lying on the specific region @varregion. The resulting object is a SpeciesOnRegion. This is useful for defining reaction schemes for MultiCompartmentReaction. This is provided to allow use of SpeciesOnRegion where Species would normally go. If the regions match, self is returned; otherwise an empty SpeciesOnRegion. """ if r == self._region(): return self else: return SpeciesOnRegion(self._species, None) @property def states(self): """A vector of all the states corresponding to this species""" all_states = node._get_states() return [all_states[i] for i in numpy.sort(self.indices())] @property def nodes(self): """A NodeList of the Node objects containing concentration data for the given Species and Region. The code node_list = ca[cyt].nodes is more efficient than the otherwise equivalent node_list = ca.nodes(cyt) because the former only creates the Node objects belonging to the restriction ca[cyt] whereas the second option constructs all Node objects belonging to the Species ca and then culls the list to only include those also belonging to the Region cyt. """ initializer._do_init() return nodelist.NodeList(itertools.chain.from_iterable([s.nodes for s in self._species()._secs if s._region == self._region()])) @property def concentration(self): """An iterable of the current concentrations.""" return self.nodes.concentration @concentration.setter def concentration(self, value): """Sets all Node objects in the restriction to the specified concentration value. This is equivalent to s.nodes.concentration = value""" self.nodes.concentration = value # 3d matrix stuff def _setup_matrices_process_neighbors(pt1, pt2, indices, euler_matrix, index, diffs, vol, areal, arear, dx): # TODO: validate this before release! is ignoring reflective boundaries the right thing to do? # (make sure that any changes here also work with boundaries that aren't really reflective, but have a 1d section attached) d = diffs[index] if pt1 in indices: ileft = indices[pt1] dleft = (d + diffs[ileft]) * 0.5 left = dleft * areal / (vol * dx) euler_matrix[index, ileft] += left euler_matrix[index, index] -= left if pt2 in indices: iright = indices[pt2] dright = (d + diffs[iright]) * 0.5 right = dright * arear / (vol * dx) euler_matrix[index, iright] += right euler_matrix[index, index] -= right # TODO: make sure that we can make this work where things diffuse across the # boundary between two regions... for the tree solver, this is complicated # because need the sections in a sorted order # ... but this is also weird syntactically, because how should it know # that two regions (e.g. apical and other_dendrites) are connected? class Species(_SpeciesMathable): def __init__(self, regions=None, d=0, name=None, charge=0, initial=None, atolscale=1): """s = rxd.Species(regions, d = 0, name = None, charge = 0, initial = None, atolscale=1) Declare a species. Parameters: regions -- a Region or list of Region objects containing the species d -- the diffusion constant of the species (optional; default is 0, i.e. non-diffusing) name -- the name of the Species; used for syncing with HOC (optional; default is none) charge -- the charge of the Species (optional; default is 0) initial -- the initial concentration or None (if None, then imports from HOC if the species is defined at finitialize, else 0) atolscale -- scale factor for absolute tolerance in variable step integrations Note: charge must match the charges specified in NMODL files for the same ion, if any. You probably want to adjust atolscale for species present at low concentrations (e.g. calcium). """ import neuron import ctypes from . import rxd # if there is a species, then rxd is being used, so we should register # this function may be safely called many times rxd._do_nbs_register() # TODO: check if "name" already defined elsewhere (if non-None) # if so, make sure other fields consistent, expand regions as appropriate self._allow_setting = True self.regions = regions self.d = d self.name = name self.charge = charge self.initial = initial self._atolscale = atolscale _all_species.append(weakref.ref(self)) # declare an update to the structure of the model (the number of differential equations has changed) neuron.nrn_dll_sym('structure_change_cnt', ctypes.c_int).value += 1 # initialize self if the rest of rxd is already initialized if initializer.is_initialized(): if _has_3d: if isinstance(regions, region.Region) and not regions._secs1d: pass elif hasattr(regions, '__len__') and all(not r._secs1d for r in regions): pass else: # TODO: remove this limitation # one strategy would be to just redo the whole thing; what are the implications of that? # (pointers would be invalid; anything else?) raise RxDException('Currently cannot add species containing 1D after 3D species defined and initialized. To work-around: reorder species definition.') self._do_init() def _do_init(self): self._do_init1() self._do_init2() self._do_init3() self._do_init4() self._do_init5() def _do_init1(self): from . import rxd # TODO: if a list of sections is passed in, make that one region # _species_count is used to create a unique _real_name for the species global _species_count regions = self._regions self._real_name = self._name initial = self.initial charge = self._charge # invalidate any old initialization of external solvers rxd._external_solver_initialized = False # TODO: check about the 0