diff --git a/doc/source/index.rst b/doc/source/index.rst index 615b23c..e535874 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -139,6 +139,15 @@ class :class:`eegmegcalc.NYHeadModel` :undoc-members: +Module :mod:`lfpykit.special` +============================= +.. automodule:: lfpykit.special + :members: + :show-inheritance: + :undoc-members: + + + Indices and tables ================== diff --git a/examples/Example_Arbor.ipynb b/examples/Example_Arbor.ipynb index 6969304..76f8a68 100644 --- a/examples/Example_Arbor.ipynb +++ b/examples/Example_Arbor.ipynb @@ -190,7 +190,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Collect CV geometry data using arbor's place_pwlin methoood\n", + "# Collect CV geometry data using arbor's place_pwlin method\n", "p = arbor.place_pwlin(morphology)\n", "x, y, z, d = [], [], [], []\n", "for m in I_m_meta:\n", @@ -354,7 +354,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.10.2" } }, "nbformat": 4, diff --git a/examples/Example_Arbor_swc.ipynb b/examples/Example_Arbor_swc.ipynb index 4e6d55c..7f28b79 100644 --- a/examples/Example_Arbor_swc.ipynb +++ b/examples/Example_Arbor_swc.ipynb @@ -10,7 +10,13 @@ "the line source approximation implementation in class `LineSourcePotential` with a passive neuron model set up in Arbor (https://arbor.readthedocs.io, https://github.com/arbor-sim/arbor). \n", "\n", "The neuron receives sinusoid synaptic current input in one arbitrary chosen control volume (CV). \n", - "Its morphology is defined in the file `single_cell.swc`" + "Its morphology is defined in the file `single_cell.swc`\n", + "\n", + "The example produces comparable output as the `Example_Neuron_swc.ipynb` and `Example_LFPy_swc.ipynb` notebooks. \n", + "\n", + "This example showcase the use of `lfpykit.special.CellGeometryArbor` and `lfpykit.special.LineSourcePotential`, \n", + "two modified classes inherited from `lfpykit.CellGeometry` and `lfpykit.LineSourcePotential`, respectively. \n", + "The modifications allow using the full geometry detail from the `.swc` morphology file for computing the extracellular potential." ] }, { @@ -225,7 +231,8 @@ "metadata": {}, "source": [ "## Compute extracellular potentials\n", - "First we define a couple of classes to interface the LFPykit library (https://LFPykit.readthedocs.io, https://github.com/LFPy/LFPykit):" + "Combining the use of `lfpykit.special.CellGeometryArbor` to represent the coordinates of segments for each CV and \n", + "`lfpykit.special.LineSourcePotential` compute the extracellular potential" ] }, { @@ -234,40 +241,8 @@ "metadata": {}, "outputs": [], "source": [ - "class ArborCellGeometry(lfpykit.CellGeometry):\n", - " '''\n", - " Class inherited from ``lfpykit.CellGeometry`` for easier forward-model\n", - " predictions in Arbor that keeps track of arbor.segment information\n", - " for each CV.\n", - "\n", - " Parameters\n", - " ----------\n", - " p: ``arbor.place_pwlin`` object\n", - " 3-d locations and cables in a morphology (cf. ``arbor.place_pwlin``)\n", - " cables: ``list``\n", - " ``list`` of corresponding ``arbor.cable`` objects where transmembrane\n", - " currents are recorded (cf. ``arbor.cable_probe_total_current_cell``)\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.CellGeometry\n", - " '''\n", - "\n", - " def __init__(self, p, cables):\n", - " x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4\n", - " c_ind = np.array([], dtype=int) # tracks which CV owns segment\n", - " for i, m in enumerate(cables):\n", - " segs = p.segments([m])\n", - " for j, seg in enumerate(segs):\n", - " x = np.row_stack([x, [seg.prox.x, seg.dist.x]])\n", - " y = np.row_stack([y, [seg.prox.y, seg.dist.y]])\n", - " z = np.row_stack([z, [seg.prox.z, seg.dist.z]])\n", - " d = np.row_stack(\n", - " [d, [seg.prox.radius * 2, seg.dist.radius * 2]])\n", - " c_ind = np.r_[c_ind, i]\n", - "\n", - " super().__init__(x=x, y=y, z=z, d=d)\n", - " self._c_ind = c_ind" + "# create ``CellGeometryArbor`` instance\n", + "cell_geometry = lfpykit.special.CellGeometryArbor(p, I_m_meta)" ] }, { @@ -276,65 +251,6 @@ "metadata": {}, "outputs": [], "source": [ - "class ArborLineSourcePotential(lfpykit.LineSourcePotential):\n", - " '''subclass of ``lfpykit.LineSourcePotential`` modified for\n", - " instances of ``ArborCellGeometry``.\n", - " Each CV may consist of several segments , and this implementation\n", - " accounts for their contributions normalized by surface area, that is,\n", - " we assume constant transmembrane current density per area across each CV\n", - " and constant current source density per unit length per segment\n", - " (inherent in the line-source approximation).\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: object\n", - " ``ArborCellGeometry`` instance or similar.\n", - " x: ndarray of floats\n", - " x-position of measurement sites (µm)\n", - " y: ndarray of floats\n", - " y-position of measurement sites (µm)\n", - " z: ndarray of floats\n", - " z-position of measurement sites (µm)\n", - " sigma: float > 0\n", - " scalar extracellular conductivity (S/m)\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.LineSourcePotential\n", - " '''\n", - "\n", - " def __init__(self, **kwargs):\n", - " super().__init__(**kwargs)\n", - " self._get_transformation_matrix = super().get_transformation_matrix\n", - "\n", - " def get_transformation_matrix(self):\n", - " '''Get linear response matrix\n", - "\n", - " Returns\n", - " -------\n", - " response_matrix: ndarray\n", - " shape (n_coords, n_cs) ndarray\n", - " '''\n", - " M_tmp = self._get_transformation_matrix()\n", - " n_cs = np.unique(self.cell._c_ind).size\n", - " M = np.zeros((self.x.size, n_cs))\n", - " for i in range(n_cs):\n", - " inds = self.cell._c_ind == i\n", - " M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] /\n", - " self.cell.area[inds].sum())\n", - "\n", - " return M" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "# create ``ArborCellGeometry`` instance\n", - "cell_geometry = ArborCellGeometry(p, I_m_meta)\n", - "\n", "# define locations where extracellular potential is predicted in vicinity\n", "# of cell.\n", "# Axis limits [x-min, x-max, y-min, y-max] (µm)\n", @@ -347,8 +263,8 @@ " int(np.diff(axis[2:]) // dz) + 1))\n", "Z = np.zeros_like(X)\n", "\n", - "# ``ArborLineSourcePotential`` instance, get mapping for all segments per CV\n", - "lsp = ArborLineSourcePotential(cell=cell_geometry,\n", + "# ``lfpykit.special.LineSourcePotential`` instance, get mapping for all segments per CV\n", + "lsp = lfpykit.special.LineSourcePotential(cell=cell_geometry,\n", " x=X.flatten(),\n", " y=Y.flatten(),\n", " z=Z.flatten())\n", @@ -368,7 +284,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -398,7 +314,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -423,7 +339,7 @@ " colors = [plt.get_cmap(cmap)(norm(v)) for v in V_m]\n", " zips = []\n", " for i in range(V_m.size):\n", - " inds = cell_geometry._c_ind == i\n", + " inds = cell_geometry._compartment_index == i\n", " zips.append(create_polygon(cell_geometry.x[inds, ].flatten(),\n", " cell_geometry.y[inds, ].flatten(),\n", " cell_geometry.d[inds, ].flatten()))\n", @@ -436,7 +352,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -464,7 +380,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ @@ -492,7 +408,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": {}, "outputs": [ { @@ -501,7 +417,7 @@ "Text(0.5, 1.0, '$V_e$ and $V_m$ at $t$=499.0 ms')" ] }, - "execution_count": 15, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, diff --git a/examples/Example_LFPy_pt3d.ipynb b/examples/Example_LFPy_pt3d.ipynb index 626bccd..dcdec23 100644 --- a/examples/Example_LFPy_pt3d.ipynb +++ b/examples/Example_LFPy_pt3d.ipynb @@ -13,7 +13,9 @@ "The example is meant to demonstrate how pt3d information (for instance the full 3D geometry information present in neuronal reconstruction file formats such as `.swc`) can be accounted for in forward-model predictions,\n", "compared to the typical LFPy setup where each compartment is treated as equivalent, straight line sources. \n", "\n", - "To deal with the pt3d information and interfacing LFPykit this notebook defines a couple inherited classes `Pt3dCellGeometry` and `Pt3dLineSourcePotential` wrapping the functionality of the corresponding LFPykit classes." + "This example showcase the use of `lfpykit.special.CellGeometryLFPyPt3d` and `lfpykit.special.LineSourcePotential`, \n", + "two modified classes inherited from `lfpykit.CellGeometry` and `lfpykit.LineSourcePotential`, respectively. \n", + "The modifications allow using the full geometry detail (pt3d points) from the morphology for computing the extracellular potential. " ] }, { @@ -39,7 +41,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -50,7 +52,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -64,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -89,118 +91,11 @@ " return cb" ] }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "class Pt3dCellGeometry(lfpykit.CellGeometry):\n", - " '''\n", - " Class inherited from ``lfpykit.CellGeometry`` for easier forward-model\n", - " predictions in LFPy that keeps track of pt3d information\n", - " for each compartment.\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: ``LFPy.Cell`` object\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.CellGeometry\n", - " '''\n", - " def __init__(self, cell):\n", - " x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4\n", - " c_ind = np.array([], dtype=int) # tracks which compartment owns segment\n", - " idx = 0\n", - " for i, sec in enumerate(cell.allseclist):\n", - " xp = (cell.arc3d[i] / cell.arc3d[i][-1])\n", - " seg_l = 1 / sec.nseg / 2 \n", - " for seg in sec:\n", - " seg_x = np.array([seg.x - seg_l, seg.x + seg_l])\n", - " ind = (xp > seg_x.min()) & (xp < seg_x.max())\n", - " seg_x = np.sort(np.r_[seg_x, xp[ind]])\n", - " \n", - " # interpolate to get values for start-, pt3d-, and end-point for search segment:\n", - " x_tmp = np.interp(seg_x, xp, cell.x3d[i])\n", - " y_tmp = np.interp(seg_x, xp, cell.y3d[i])\n", - " z_tmp = np.interp(seg_x, xp, cell.z3d[i])\n", - " d_tmp = np.interp(seg_x, xp, cell.diam3d[i])\n", - " \n", - " # store\n", - " x = np.row_stack([x, np.c_[x_tmp[:-1], x_tmp[1:]]])\n", - " y = np.row_stack([y, np.c_[y_tmp[:-1], y_tmp[1:]]])\n", - " z = np.row_stack([z, np.c_[z_tmp[:-1], z_tmp[1:]]])\n", - " d = np.row_stack([d, np.c_[d_tmp[:-1], d_tmp[1:]]])\n", - " \n", - " c_ind = np.r_[c_ind, np.ones(seg_x.size - 1) * idx]\n", - " idx += 1\n", - "\n", - " super().__init__(x=x, y=y, z=z, d=d)\n", - " self._c_ind = c_ind" - ] - }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], - "source": [ - "class Pt3dLineSourcePotential(lfpykit.LineSourcePotential):\n", - " '''subclass of ``lfpykit.LineSourcePotential`` modified for\n", - " instances of ``Pt3dCellGeometry``.\n", - " Each compartment may consist of several segments, and this implementation\n", - " accounts for their contributions normalized by surface area, that is,\n", - " we assume constant transmembrane current density per area across each compartment\n", - " and constant current source density per unit length per segment\n", - " (inherent in the line-source approximation).\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: object\n", - " ``Pt3dCellGeometry`` instance or similar.\n", - " x: ndarray of floats\n", - " x-position of measurement sites (µm)\n", - " y: ndarray of floats\n", - " y-position of measurement sites (µm)\n", - " z: ndarray of floats\n", - " z-position of measurement sites (µm)\n", - " sigma: float > 0\n", - " scalar extracellular conductivity (S/m)\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.LineSourcePotential\n", - " '''\n", - "\n", - " def __init__(self, **kwargs):\n", - " super().__init__(**kwargs)\n", - " self._get_transformation_matrix = super().get_transformation_matrix\n", - "\n", - " def get_transformation_matrix(self):\n", - " '''Get linear response matrix\n", - "\n", - " Returns\n", - " -------\n", - " response_matrix: ndarray\n", - " shape (n_coords, n_CVs) ndarray\n", - " '''\n", - " M_tmp = self._get_transformation_matrix()\n", - " n_comps = np.unique(self.cell._c_ind).size\n", - " M = np.zeros((self.x.size, n_comps))\n", - " for i in range(n_comps):\n", - " inds = self.cell._c_ind == i\n", - " M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] /\n", - " self.cell.area[inds].sum())\n", - "\n", - " return M" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], "source": [ "# set seed\n", "np.random.seed(1234)" @@ -208,7 +103,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -217,7 +112,7 @@ "1.0" ] }, - "execution_count": 8, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -259,7 +154,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -290,7 +185,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -298,9 +193,9 @@ "# run simulation, compute measurement\n", "###############################################################################\n", "cell.simulate(variable_dt=True,\n", - " rec_imem=True,\n", - " rec_vmem=True,\n", - " )\n", + " rec_imem=True,\n", + " rec_vmem=True,\n", + " )\n", "\n", "# copy some output to GeometryCell for plotting\n", "cell.time = cell.tvec\n", @@ -310,7 +205,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -336,7 +231,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -344,10 +239,10 @@ "# compute extracellular potential using pt3d information\n", "###############################################################################\n", "# create ``Pt3dCellGeometry`` instance\n", - "cell_geometry = Pt3dCellGeometry(cell)\n", + "cell_geometry = lfpykit.special.CellGeometryLFPyPt3d(cell)\n", "\n", "# LineSourcePotential object, get mapping\n", - "lsp_pt3d = Pt3dLineSourcePotential(cell=cell_geometry, x=x, y=y, z=z)\n", + "lsp_pt3d = lfpykit.special.LineSourcePotential(cell=cell_geometry, x=x, y=y, z=z)\n", "M_pt3d = lsp_pt3d.get_transformation_matrix()\n", "\n", "# Extracellular potential using pt3d information at last time step \n", @@ -357,7 +252,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [ { diff --git a/examples/Example_LFPy_swc.ipynb b/examples/Example_LFPy_swc.ipynb index a7e579b..918325a 100644 --- a/examples/Example_LFPy_swc.ipynb +++ b/examples/Example_LFPy_swc.ipynb @@ -12,8 +12,11 @@ "The neuron receives sinusoid synaptic current input in one arbitrary chosen segment. \n", "Its morphology is defined in the file `single_cell.swc`. \n", "\n", - "The example produces comparable output as the `Example_Arbor_swc.ipynb` notebook. \n", - "For this a couple of modified `CellGeometry` and `LineSourcePotential` classes have been defined below, similar to `Example_LFPy_pt3d.ipynb`" + "The example produces comparable output as the `Example_Arbor_swc.ipynb` and `Example_Neuron_swc.ipynb` notebooks. \n", + "\n", + "This example showcase the use of `lfpykit.special.CellGeometryLFPyPt3d` and `lfpykit.special.LineSourcePotential`, \n", + "two modified classes inherited from `lfpykit.CellGeometry` and `lfpykit.LineSourcePotential`, respectively. \n", + "The modifications allow using the full geometry detail from the `.swc` morphology file for computing the extracellular potential. " ] }, { @@ -67,113 +70,6 @@ "execution_count": 5, "metadata": {}, "outputs": [], - "source": [ - "class Pt3dCellGeometry(lfpykit.CellGeometry):\n", - " '''\n", - " Class inherited from ``lfpykit.CellGeometry`` for easier forward-model\n", - " predictions in LFPy that keeps track of pt3d information\n", - " for each compartment.\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: ``LFPy.Cell`` object\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.CellGeometry\n", - " '''\n", - " def __init__(self, cell):\n", - " x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4\n", - " c_ind = np.array([], dtype=int) # tracks which compartment owns segment\n", - " idx = 0\n", - " for i, sec in enumerate(cell.allseclist):\n", - " xp = (cell.arc3d[i] / cell.arc3d[i][-1])\n", - " seg_l = 1 / sec.nseg / 2 \n", - " for seg in sec:\n", - " seg_x = np.array([seg.x - seg_l, seg.x + seg_l])\n", - " ind = (xp > seg_x.min()) & (xp < seg_x.max())\n", - " seg_x = np.sort(np.r_[seg_x, xp[ind]])\n", - " \n", - " # interpolate to get values for start-, pt3d-, and end-point for search segment:\n", - " x_tmp = np.interp(seg_x, xp, cell.x3d[i])\n", - " y_tmp = np.interp(seg_x, xp, cell.y3d[i])\n", - " z_tmp = np.interp(seg_x, xp, cell.z3d[i])\n", - " d_tmp = np.interp(seg_x, xp, cell.diam3d[i])\n", - " \n", - " # store\n", - " x = np.row_stack([x, np.c_[x_tmp[:-1], x_tmp[1:]]])\n", - " y = np.row_stack([y, np.c_[y_tmp[:-1], y_tmp[1:]]])\n", - " z = np.row_stack([z, np.c_[z_tmp[:-1], z_tmp[1:]]])\n", - " d = np.row_stack([d, np.c_[d_tmp[:-1], d_tmp[1:]]])\n", - " \n", - " c_ind = np.r_[c_ind, np.ones(seg_x.size - 1) * idx]\n", - " idx += 1\n", - "\n", - " super().__init__(x=x, y=y, z=z, d=d)\n", - " self._c_ind = c_ind" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "class Pt3dLineSourcePotential(lfpykit.LineSourcePotential):\n", - " '''subclass of ``lfpykit.LineSourcePotential`` modified for\n", - " instances of ``Pt3dCellGeometry``.\n", - " Each compartment may consist of several segments, and this implementation\n", - " accounts for their contributions normalized by surface area, that is,\n", - " we assume constant transmembrane current density per area across each compartment\n", - " and constant current source density per unit length per segment\n", - " (inherent in the line-source approximation).\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: object\n", - " ``Pt3dCellGeometry`` instance or similar.\n", - " x: ndarray of floats\n", - " x-position of measurement sites (µm)\n", - " y: ndarray of floats\n", - " y-position of measurement sites (µm)\n", - " z: ndarray of floats\n", - " z-position of measurement sites (µm)\n", - " sigma: float > 0\n", - " scalar extracellular conductivity (S/m)\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.LineSourcePotential\n", - " '''\n", - "\n", - " def __init__(self, **kwargs):\n", - " super().__init__(**kwargs)\n", - " self._get_transformation_matrix = super().get_transformation_matrix\n", - "\n", - " def get_transformation_matrix(self):\n", - " '''Get linear response matrix\n", - "\n", - " Returns\n", - " -------\n", - " response_matrix: ndarray\n", - " shape (n_coords, n_CVs) ndarray\n", - " '''\n", - " M_tmp = self._get_transformation_matrix()\n", - " n_comps = np.unique(self.cell._c_ind).size\n", - " M = np.zeros((self.x.size, n_comps))\n", - " for i in range(n_comps):\n", - " inds = self.cell._c_ind == i\n", - " M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] /\n", - " self.cell.area[inds].sum())\n", - "\n", - " return M" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], "source": [ "def custom_fun(cell, nseg=3):\n", " # approximate the same cable discretization as the Arbor example w. nseg=3 per branch\n", @@ -187,7 +83,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -221,16 +117,31 @@ "cell.simulate(rec_imem=True, rec_vmem=True)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## compute extracellular potentials\n", + "Combining the use of `lfpykit.special.CellGeometryArbor` to represent the coordinates of `pt3d`-points for each compartment and \n", + "`lfpykit.special.LineSourcePotential` compute the extracellular potential." + ] + }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# create ``lfpykit.special.CellGeometryLFPyPt3d`` instance\n", + "cell_geometry = lfpykit.special.CellGeometryLFPyPt3d(cell)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ - "###############################################################################\n", - "# compute extracellular potential using pt3d information\n", - "###############################################################################\n", - "\n", "# locations where extracellular potential is predicted \n", "dx = 1\n", "dz = 1\n", @@ -240,77 +151,18 @@ " np.linspace(axis[2], axis[3], int(np.diff(axis[2:]) // dz) + 1))\n", "Z = np.zeros_like(X)\n", "\n", - "# create ``Pt3dCellGeometry`` instance\n", - "cell_geometry = Pt3dCellGeometry(cell)\n", "\n", "# LineSourcePotential object, get mapping\n", - "lsp_pt3d = Pt3dLineSourcePotential(cell=cell_geometry, \n", - " x=X.flatten(), \n", - " y=Y.flatten(), \n", - " z=Z.flatten())\n", - "M = lsp_pt3d.get_transformation_matrix()\n", + "lsp = lfpykit.special.LineSourcePotential(cell=cell_geometry, \n", + " x=X.flatten(), \n", + " y=Y.flatten(), \n", + " z=Z.flatten())\n", + "M = lsp.get_transformation_matrix()\n", "\n", "# Extracellular potential in x,z-plane coordinates \n", "V_e = M @ cell.imem" ] }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(-143.00076696500028,\n", - " 363.01610626500576,\n", - " -76.67346928554844,\n", - " 66.70703218377005)" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAABIwAAAHSCAYAAACO6lCPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJvklEQVR4nO3deVDXd57v+9eHHRFEBRc2AUVkERdwiRsatywdE5OYxCQmMS6d7jh1761zb812686ZOXXrzpxzZqbOuTPdc6N20k62TqfN0nZPC7iwiEtwQwVBxQXFDQURkfX3uX+E0P7SGDd+fFmejyoq8P2xvLHqK/jM9/v+GWutAAAAAAAAgO94OT0AAAAAAAAAehaCEQAAAAAAANwQjAAAAAAAAOCGYAQAAAAAAAA3BCMAAAAAAAC4IRgBAAAAAADAjY/TA9yvsLAwGxsb6/QYAAAAAAAAfcb+/furrbXh3z/ea4JRbGysioqKnB4DAAAAAACgzzDGnO3sOLekAQAAAAAAwA3BCAAAAAAAAG4IRgAAAAAAAHBDMAIAAAAAAIAbghEAAAAAAADcEIwAAAAAAADghmAEAAAAAAAANwQjAAAAAAAAuCEYAQAAAAAAwA3BCAAAAAAAAG4IRgAAAAAAAHBDMAIAAAAAAIAbghEAAAAAAADcEIwAAAAAAADghmAEAAAAAAAANwQjAAAAAAAAuCEYAQAAAAAAwA3BCAAAAAAAAG4IRgAAAAAAAHBDMAIAAAAAAIAbghEAAAAAAADcEIwAAAAAAADghmAEAAAAAAAANwQjAAAAAAAAuCEYAQAAAAAAwA3BCAAAAAAAAG4IRgAAAAAAAHBDMAIAAAAAAIAbghEAAAAAAADcEIwAAAAAAADghmAEAAAAAAAANwQjAAAAAAAAuCEYAQAAAAAAwA3BCAAAAAAAAG4IRgAAAAAAAHDj48lPboxJlPSrOw7FS/q/JIVKWiPpavvxv7LW/t6TswAAAAAAAOD+eDQYWWvLJE2UJGOMt6QLkr6QtFLSP1tr/7snvz4AAAAAAAAeXHfekjZf0ilr7dlu/JoAAAAAAAB4QN0ZjF6R9Mkdb68zxhQbY35hjBncjXMAAAAAAADgB3RLMDLG+ElaIunX7Yd+Lmm0vr1d7aKkf7zLx601xhQZY4quXr3a2bsAAAAAAACgi3XXFUZPSjpgrb0sSdbay9baNmutS9J6SVM7+yBr7XvW2gxrbUZ4eHg3jQoAAAAAANC/dVcwWq47bkczxoy847Glko520xwAAAAAAAC4B48+S5okGWMGSFoo6cd3HP6vxpiJkqykM997DAAAAAAAAA7yeDCy1jZIGvq9Yys8/XUBAAAAAADwcLrzWdIAAAAAAADQCxCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADc+Tg8AAAAAAHhwLS0tOnHihEpKSuTr66tnn33W6ZEA9CEEIwAAAADoJZqbmzsiUUlJScfxAQMGyForY4yD0wHoSwhGAAAAANCDNTU1qby8XCUlJTp+/LjbY8nJyUpOTlZCQgKxCECXIhgBAAAAQA/T2NjYEYnKysrcHhs3bpxSUlI0duxY+fn5OTQhgL6OYAQAAAAAPUBjY6OOHz+u0tJSlZeXuz2WmJio5ORkJSYmyt/f36EJAfQnBCMAAAAAcMjt27d1/PhxlZSU6NSpU7LWdjw2duzYjkgUEBDg4JQA+iOCEQAAAAB0o4aGho5IdPr0ablcro7HxowZo5SUFCUmJiowMNDBKQH0dwQjAAAAAPCwW7duqbS0VKWlpTpz5oystR3PajZ69GglJycrKSmJSASgxyAYAQAAAIAH1NfXq7S0VCUlJaqsrJQxpuOWs9jYWKWkpCgpKUkDBgxweFIA+FMEIwAAAADoIjdv3uyIRBcuXJCXl1dHJIqKiuqIREFBQQ5PCgA/jGAEAAAAAI+grq5OJSUlKi0t1cWLF+Xj49MRiUaOHKnk5GQlJydr4MCBDk8KAPePYAQAAAAAD+jGjRsqKSlRSUmJrly5Ij8/v46dROHh4R2RKDg42OlRAeChEIwAAAAA4D7U1tZ2RKLq6mr5+/t3XEkUGhqqlJQUJScnKyQkxOFJAeDREYwAAAAA4C5qamo6ItH169cVGBjYcSVRcHBwRyQaNGiQ06MCQJciGAEAAADAHa5du9YRierq6jRgwAC5XC4ZYxQYGKjk5GSlpKQoNDTU6VEBwGMIRgAAAAD6verq6o5IVF9fr6CgoI7bzXx8fDquJBoyZIjDkwJA9yAYAQAAAOiXrl69qmPHjqm0tFQNDQ0KDg6WMUaSZIxRSkqKUlJSNHToUIcnBYDuRzACAAAA0C9Ya3XlypWOK4mampoUEhIiX19fGWPU1tbWcbtZWFiY0+MCgKMIRgAAAAD6LGutLl++3BGJWlpaFBoaqsDAQDU3N6upqanjdrNhw4Y5PS4A9BgeD0bGmDOSbkpqk9Rqrc0wxgyR9CtJsZLOSHrJWlvj6VkAAAAA9H3WWl26dKnjdjOXy6XBgwcrODhYNTU1qq+vV0pKilJTUxUeHt5xGxoA4I+66wqjedba6jve/gtJ26y1f2+M+Yv2t/+8m2YBAAAA0MdYa1VVVaWSkhKVlpZKkoYMGaLBgwfr+vXrqq2tVUpKihYvXqzhw4cTiQDgHpy6Je1ZSXPbX/+lpJ0iGAEAAAB4ANZaXbhwoeN2M29vb4WFhSksLEzXrl3TtWvXlJKSovnz52vEiBFEIgB4AN0RjKykLGOMlfT/WWvfkzTcWntRkqy1F40x3CwMAAAA4J6stTp//nzH7WZ+fn4KDw/XyJEjVV1drcuXLys5OVmZmZmKiIggEgHAQ+qOYDTTWlvVHoWyjTHH7/cDjTFrJa2VpJiYGE/NBwAAAKAHs9bq3LlzHbebBQQEaNiwYYqOjtaVK1d04cIFJScna+bMmYqMjCQSAUAX8HgwstZWtf/3ijHmC0lTJV02xoxsv7popKQrd/nY9yS9J0kZGRnW07MCAAAA6BlcLpdbJBowYIBGjBih2NhYXb58WefOnVNycrKmTZumqKgoIhEAdDGPBiNjTJAkL2vtzfbXF0n6O0lfS3pT0t+3//crT84BAAAAoOdzuVw6e/ZsRyQKDg7WiBEjNHr0aF2+fFkVFRVKSkrSU089pZiYGCIRAHiQp68wGi7pi/a/yH0kfWyt/YMx5htJnxljVkk6J2mZh+cAAAAA0AO5XC6dPn1aJSUlOn78uAYNGqSRI0cqMTFRFy9e1IkTJ5SUlKTFixcrJiZGXl5eTo8MAP2CsbZ33OmVkZFhi4qKnB4DAAAAwCNqa2vriERlZWUKDQ1VZGSkrLWqqqpSbW2tkpKSlJKSolGjRhGJAMCDjDH7rbUZ3z/eHUuvAQAAAPRzbW1tqqio6IhEQ4cOVWRkpMaPH68LFy7o6NGjGjdunObPn6/Y2FgiEQA4jGAEAAAAwCNaW1t16tQplZaWqqysTOHh4YqOjtbEiRN1/vx5FRcXKzExUZmZmYqLi5O3t7fTIwMA2hGMAAAAAHSZlpYWnTp1SiUlJTpx4oSGDx+u6Ohopaenq7KyUgcOHFBiYqJmz56t+Ph4IhEA9FAEIwAAAACPpKWlRSdPnuyIRCNHjlRMTIymTp2qc+fOqaioSGPHjtXMmTMVHx8vHx/+GQIAPR1/UwMAAAB4YM3NzTpx4oRKS0t18uRJRUREKDY2VmFhYTp79qz27t2rsWPHatq0aRozZgyRCAB6Gf7WBgAAAHBfmpubVV5erpKSElVUVCgqKkqxsbEaPny4zpw5o8LCQo0ZM0ZTpkzRmDFj5Ovr6/TIAICHRDACAAAAcFdNTU1ukSgmJkZxcXGKiorS6dOnVVBQoNGjR2vy5Ml65ZVXiEQA0EcQjAAAAAC4aWxsVFlZmUpLS3X69GmNGjVK8fHxGjVqlCoqKpSXl6f4+HhNmDBBy5Ytk5+fn9MjAwC6GMEIAAAAgG7fvq2ysjKVlJTo7Nmzio2N1ejRoxUfH69Tp05px44diouL0/jx4/XCCy/I39/f6ZEBAB5EMAIAAAD6qYaGho5IdO7cOcXHx2vs2LEaO3asTp06pe3bt2vUqFFKTk7W0qVLFRAQ4PTIAIBuQjACAAAA+pGGhgaVlpaqtLRU58+fV3x8vJKSkpSUlKSTJ08qJydH0dHRSklJ0bPPPkskAoB+imAEAAAA9HG3bt1SaWmpSkpKVFVVpdGjRys1NVWpqak6ceKEsrKyFBUVpZSUFD3zzDMKDAx0emQAgMMIRgAAAEAfVF9f3xGJLl68qISEBE2cOFETJ07UiRMn9Ic//EGRkZFKTk7W008/rQEDBjg9MgCgByEYAQAAAH3EzZs3VVJSotLSUl2+fFkJCQlKT0+XJJWXl+v3v/+9IiIilJycrCeeeEJBQUEOTwwA6KkIRgAAAEAvVldXp5KSEpWUlOjq1asaO3aspk6dKkkqKyvT7373Ow0fPlwpKSlatGiRBg4c6PDEAIDegGAEAAAA9DJ1dXU6duyYSkpKdO3aNSUmJmrGjBmSvo1Ev/3tbzVs2DClpKRowYIFCg4OdnhiAEBvQzACAAAAeomGhgbl5eWpuLhYiYmJmjVrlowxKi0t1ddff62wsDAlJydr3rx5CgkJcXpcAEAvRjACAAAAerjW1lbt27dPu3btUnJysl599VUdOnRIX3/9tQYPHqyUlBTNnTtXgwYNcnpUAEAfQTACAAAAeihrrUpKSpSTk6Nhw4bp1VdfVWlpqT7++GNNmTJFa9asUWhoqNNjAgD6IIIRAAAA0ANVVlYqKytLra2t+tGPfqQbN27o008/1ejRo/XOO+9wyxkAwKMIRgAAAEAPUlNTo5ycHJ0/f75jF1FWVpb8/f21fPlyRUREOD0iAKAfIBgBAAAAPcDt27eVl5enw4cPa/r06ZozZ4527Nihy5cva+HChUpKSpIxxukxAQD9BMEIAAAAcFBbW5u++eYb5efna9y4cXr77be1f/9+/fKXv9SMGTP04osvyseHX9sBAN2LnzwAAACAA6y1Ki0tVU5OjoYOHaoVK1aosrJSH3zwgRITE/XTn/5UAwcOdHpMAEA/RTACAAAAutn58+eVlZWl5uZmPf3007LWavPmzQoKCtLrr7+uESNGOD0iAKCfIxgBAAAA3aSmpkbbt2/X2bNnNW/ePEVGRionJ0fXrl3TwoULlZiYyJ4iAECPQDACAAAAPKyxsVH5+fk6ePCgpk6dqgULFqiwsFA5OTmaNWuWXn75ZXl7ezs9JgAAHbycHgAAAADoq9ra2rR37179y7/8i27fvq0f//jHCggI0HvvvSeXy6V3331Xjz32GLEIANDjcIURAAAA0MWstSorK1N2drYGDx6sFStW6MaNG/r3f/93hYaG6s0339SwYcOcHhMAgLsiGAEAAABdqKqqSllZWWpoaNCTTz6p4OBgZWVlqa6uTosXL1ZCQoLTIwIAcE8EIwAAAKAL1NbWavv27Tp9+rTmzZunhIQE5ebmqrS0VJmZmUpPT+fWMwBAr0EwAgAAAB5BY2OjCgoKdODAAU2ZMkVPPPGEDh48qJ///OdKS0vTunXrFBgY6PSYAAA8EIIRAAAA8BDa2tp04MAB5ebmKiEhQe+8844uXLigDRs2KDw8XG+//bbCwsKcHhMAgIdCMAIAAAAegLVW5eXlys7OVkhIiF5//XVZa7V582bdvn1bP/rRjxQfH+/0mAAAPBKCEQAAAHCfLl68qKysLNXX12vx4sUaMWKEtm/frhMnTmju3LmaPHmyvLy8nB4TAIBHRjACAAAA7uHGjRvavn27KioqlJmZqbS0NO3du1dffPGFJk2apHXr1ikgIMDpMQEA6DIEIwAAAOAumpqaVFBQoP379ysjI0Pr1q3TiRMn9LOf/UwRERFavXq1hgwZ4vSYAAB0OYIRAAAA8D0ul6tjofXo0aP14x//WPX19froo4/U0tKi5557TrGxsU6PCQCAxxCMAAAAgHbWWp04cULZ2dkaOHCgXn31VQUFBWnbtm06ffq05s2bpwkTJrCnCADQ5xGMAAAAAEmXLl1SVlaW6urqtHDhQsXFxamwsFD79u3ruB3Nz8/P6TEBAOgWBCMAAAD0a3V1ddqxY4dOnDihzMxMTZ48WUePHtW//uu/KiYmRmvXrlVoaKjTYwIA0K0IRgAAAOiXmpubtWvXLn3zzTeaPHmy1q1bpytXrugXv/iFjDF68cUXFR0d7fSYAAA4gmAEAACAfsXlcungwYPauXOn4uLitHbtWknSli1bVFlZqfnz52v8+PEyxjg8KQAAziEYAQAAoN84efKksrOzFRgYqOXLl2vo0KEqKCjQ/v37NW3aND377LPy9fV1ekwAABxHMAIAAECfd/nyZWVnZ6umpkYLFy7U2LFjdejQIX3yyScaPXq03nnnHYWEhDg9JgAAPQbBCAAAAH3WzZs3tWPHDpWXl2v27NnKyMhQZWWl1q9fL19fX73yyiuKjIx0ekwAAHocghEAAAD6nObmZhUWFmrfvn2aNGmS1q1bp4aGBn3++ee6dOmSFixYoOTkZPYUAQBwFx4NRsaYaEmbJI2Q5JL0nrX2fxhj/rOkNZKutr/rX1lrf+/JWQAAAND3uVwuHT58WDt27NCoUaO0du1aBQQEKC8vT4cOHdKMGTP0wgsvyMeH/28KAMAP8fRPylZJ/8lae8AYEyxpvzEmu/2xf7bW/ncPf30AAAD0E6dOnVJ2drb8/Pz00ksvKSIiQvv371dubq7Gjh2rn/70pxo4cKDTYwIA0Ct4NBhZay9Kutj++k1jTKkkbhIHAABAl7ly5Yqys7N17do1LVy4UOPGjVNFRYX+7d/+TUFBQXr99dc1YsQIp8cEAKBX6bZrcY0xsZImSdoraaakdcaYNyQV6durkGo6+Zi1ktZKUkxMTHeNCgAAgF6gvr5eO3bs0PHjxzV79my98sorqqmp0SeffKLq6motWrRIiYmJ7CkCAOAhGGut57+IMQMl5Ur6v621m40xwyVVS7KS/oukkdbat3/oc2RkZNiioiKPzwoAAICeraWlRbt379aePXs0YcIEzZkzR5K0c+dOHT16VDNnztTUqVPZUwQAwH0wxuy31mZ8/7jHf4oaY3wl/UbSR9bazZJkrb18x+PrJW3x9BwAAADo3ay1HQuto6OjtWbNGoWEhKioqEh5eXlKTk7WT3/6UwUFBTk9KgAAvZ6nnyXNSNooqdRa+093HB/Zvt9IkpZKOurJOQAAANC7nT59WllZWfLx8dGLL76oqKgonThxQh999JFCQ0P15ptvatiwYU6PCQBAn+HpK4xmSloh6Ygx5lD7sb+StNwYM1Hf3pJ2RtKPPTwHAAAAeqGrV68qJydHV69e1fz585WcnKyrV6/qww8/1I0bN7R48WKNGTOGPUUAAHQxTz9LWoGkzn56/96TXxcAAAC9261bt7Rz506VlJRo1qxZWrZsmZqamvS73/1OpaWlmjNnjjIyMuTt7e30qAAA9ElsAgQAAECP0dLSoj179mj37t1KS0vTu+++Kz8/P+3bt0+7du3S+PHjtW7dOgUGBjo9KgAAfRrBCAAAAI6z1qq4uFjbt29XZGSkVq9ercGDB+v48ePKzs5WeHi4Vq5cqbCwMKdHBQCgXyAYAQAAwFFnzpxRVlaWvLy89MILLygmJkaXLl3Spk2b1NDQoKefflqjR492ekwAAPoVghEAAAAcUV1drZycHF2+fFnz589XSkqKbt26pa+//lrl5eWaO3euJk+eLC8vL6dHBQCg3yEYAQAAoFvdunVLubm5Onr0qGbOnKkXX3xRklRQUKDdu3dr0qRJWrdunQICAhyeFACA/otgBAAAgG7R2tqqPXv2qLCw0G159bFjx5STk6ORI0dq9erVGjJkiNOjAgDQ7xGMAAAA4FHWWh09elTbtm3TyJEjtWrVKg0dOlQXLlzQp59+qpaWFj333HOKjY11elQAANCOYAQAAACPOXv2rLKysiSpIwrV1dXpiy++UEVFhR5//HFNmDCBPUUAAPQwBCMAAAB0uWvXriknJ0cXL17U/PnzlZqaqtbWVu3cuVP79u1Tenq61q1bJ39/f6dHBQAAnSAYAQAAoMs0NDQoNzdXR44c0YwZM/T888/Lx8dHR44c0bZt2xQdHa21a9cqNDTU6VEBAMAPIBgBAADgkbW2tmrfvn3atWuXkpOT9e677yooKEiVlZX6wx/+IEl64YUXFBMT4/CkAADgfhCMAAAA8NCstTp27Ji2bdumYcOGaeXKlQoLC1Ntba0+//xzVVZWav78+Ro/fryMMU6PCwAA7hPBCAAAAA+lsrJSW7dulcvl0pIlSxQXF6empiZt27ZN+/fv19SpU7VkyRL5+fk5PSoAAHhABCMAAAA8kOvXr2vbtm06f/68Hn/8caWlpclaq4MHD2rHjh2Ki4vTO++8o5CQEKdHBQAAD4lgBAAAgPty+/Zt5eXl6fDhw5o+fbqee+45+fr66syZM9q6dat8fX318ssvKzIy0ulRAQDAIyIYAQAA4Ae1trbqm2++UUFBgZKSkvTTn/5UAwcO1PXr15WTk6OqqiotWLBAKSkp7CkCAKCPIBgBAACgU9ZalZaWKicnR2FhYXrrrbcUHh6uxsZGZWdn6+DBg3rssce0dOlS+fr6Oj0uAADoQgQjAAAA/Inz588rKytLzc3N+tGPfqT4+Hi5XC4VFRVp586dSkhI0E9+8hMFBwc7PSoAAPAAghEAAAA61NTUaNu2bTp37lzHQmsvLy+dOnVKWVlZCgwM1GuvvaaRI0c6PSoAAPAgghEAAAB0+/Zt5efn69ChQ5o2bZqWLFkiPz8/VVdXKzs7W1evXtXChQs1btw49hQBANAPEIwAAAD6sba2NhUVFSk/P1+JiYkdt5ndvn1bf/jDH1RcXKxZs2Zp2bJl8vHhV0cAAPoLfuoDAAD0Q9ZaHT9+XDk5ORoyZIhWrFih4cOHq62tTXv37lVeXp6SkpL07rvvKigoyOlxAQBANyMYAQAA9DMXLlxQVlaWGhsb9dRTT2n06NGy1urEiRPKyspSSEiI3njjDQ0fPtzpUQEAgEMIRgAAAP1EbW2ttm/frjNnzmju3LmaOHGivLy8dOXKFWVlZam2tlaLFi1SQkICe4oAAOjnCEYAAAB9XGNjo/Lz83Xw4EFNnTpVP/rRj+Tn56dbt25p586dKikp0ezZszVlyhR5e3s7PS4AAOgBCEYAAAB9VFtbm/bv36+8vDwlJCR0LLRua2vT7t27VVBQoNTUVK1bt06BgYFOjwsAAHoQghEAAEAfY61VWVmZcnJyNGjQIL3++usaMWJEx6Lr7OxsDR06VCtXrlRYWJjT4wIAgB6IYAQAANCHVFVVKSsrSw0NDVq8eLHGjBkjY4wuXbqkrVu36tatWx2LrgEAAO6GYAQAANAH3LhxQ9u3b1dFRYXmzp2rSZMmycvLS/X19dq+fbvKy8uVmZmp9PR0eXl5OT0uAADo4QhGAAAAvVhTU5MKCgq0f/9+ZWRkaN26dfL391dra6sKCwtVWFioiRMnat26dQoICHB6XAAA0EsQjAAAAHohl8ul/fv3Kzc3V2PGjNE777yjkJAQWWt17Ngx5eTkaMSIEVq9erWGDBni9LgAAKCXIRgBAAD0ItZanThxQtnZ2QoODtZrr72mkSNHSvp2f9HWrVvV1NSkJUuWKC4uzuFpAQBAb0UwAgAA6CUuXryo7Oxs3bx5UwsXLlRCQoKMMaqrq9P27dt16tQpzZs3TxMnTmRPEQAAeCQEIwAAgB7uziCUmZmpyZMny8vLSy0tLSosLNTevXuVnp7esb8IAADgURGMAAAAeqimpibt2rVLRUVFbkHIWqvi4mJt27ZN0dHRWrt2rUJDQ50eFwAA9CEEIwAAgB7G5XLp4MGD2rlzp+Lj4/XjH/9YgwYNkiRVVlZq69atstbqhRdeUExMjMPTAgCAvohgBAAA0ENYa3Xy5EllZ2drwIABWr58uSIiIiRJtbW12rZtm86ePav58+crLS1NxhiHJwYAAH0VwQgAAKAHuHTpkrKzs3Xjxg0tXLhQY8eOlTFGzc3NKigoUFFRkaZMmaJnnnlGfn5+To8LAAD6OIIRAACAg27evKnt27frxIkTmjNnjtLT0+Xt7S1rrQ4ePKgdO3YoNjbW7bY0AAAATyMYAQAAOKC5uVmFhYXat2+fJk+erHXr1ikgIECSdPbsWW3dulXe3t566aWXFBUV5fC0AACgvyEYAQAAdCOXy6VDhw5p586dio2NdXuGs5qaGmVnZ6uqqkoLFixQSkoKe4oAAIAjCEYAAADd5NSpU8rKylJAQIBefvllRUZGSpKampqUl5engwcPavr06Vq6dKl8fX0dnhYAAPRnBCMAAAAPu3LlirKyslRTU6MFCxZo3LhxMsbI5XLp4MGD2rlzp8aMGaOf/OQnCg4OdnpcAAAAghEAAICn1NfXa/v27SorK9OcOXOUkZEhb29vSVJFRYW2bt2qwMBAvfrqqxo5cqTD0wIAAPwRwQgAAKCLNTc3a/fu3dq7d68mTpyoP/uzP+tYaH3t2jVlZWXp6tWrWrhwYcfVRgAAAD0JwQhAn/YP//APamxsdHqMPisgIEB//ud/7vQYQI/hcrlUXFys7du3KyYmRmvWrNHgwYMlSbdv31Zubq6Ki4s1c+ZMLVu2TD4+/CoGAAB6Jsd+SzHGPCHpf0jylrTBWvv3Ts0CoO9qbGzU3/zN3zg9Rp/1t3/7t06PAPQYFRUVysrKkp+fn1566SVFRUVJktra2rR//37l5eVp3LhxevfddxUUFOTwtAAAAD/MkWBkjPGW9K+SFko6L+kbY8zX1toSJ+YBADy4CxcuSJL+6Z/+yeFJ+i5jjEJCQhQaGvonL4MGDeLqlB7iypUrysnJUXV1tRYsWKCkpKSOW8xOnDihrKwshYSEaMWKFRo+fLjD0wIAANwfp37TnCrppLW2QpKMMZ9KelZSnw9G3B4DoK+4efOmJGn16tUOT9J3uVwu1dXVqba2VrW1tbpw4YKOHTum2tpa1dXVKTAwsNOYRFDqHvX19dq5c6dKS0s1a9YsvfTSSx1/5levXu14VrRFixYpISGBPUUAAKBXceo3yUhJlXe8fV7StO+/kzFmraS1khQTE9M9k3kYt8cA3YtbpjwvJCTE6RH6tNDQ0E5/BrpcLtXX13fEJIJS92lpadHu3bu1Z88eTZgwQevWrVNgYKAkqaGhQTt37tSxY8c0e/ZsTZkypeNZ0QAAAHoTp35T7Ox/sdk/OWDte5Lek6SMjIw/ebw3CggI4B+wAPqUQ4cOOT1Cn+Xj46NBgwYpNDRUAwcOdLtCxcvLSyEhIQoJCemSoPTd17nzhaDkzlrbsdA6KipKq1ev1pAhQyR9u6do3759KigoUGpqqt59910NGDDA4YkBAAAenlO/CZ6XFH3H21GSqhyapVvxbEJA9yLQes53u1jOnDnj7CB9WEtLi27cuKHa2lo1NjZ2GnW+e+mKoFRVVaWSkhKCUidOnz6trKws+fj46MUXX1R09Le/xlhrVVZWpuzsbA0dOlRvvfWWwsPDHZ4WAADg0Tn1m943khKMMXGSLkh6RdKrDs0CAHgI3z1V+HPPPefsIP3EnfHou5eysrKO13taUBo0aJB8fX2784/II6qrq5Wdna0rV65o/vz5SklJ6fhzvHz5srZu3ar6+no9+eSTGjNmjMPTAgAAdB1HgpG1ttUYs07SVknekn5hrT3mxCwAAPQGvr6+CgsLU1hYWKePE5S61q1bt7Rz506VlJRo5syZWrZsWccVVfX19dqxY4fKysqUmZmp9PR0eXl5OTwxAABA1zLW9o7VQBkZGbaoqMjpMQD0MjwzoWcFBARwq20v0VlQuvPlQYPSvXQWlO586alBqaWlRXv37lVhYaHS0tI0Z86cjl1Era2t2rNnjwoLCzVhwgRlZmYqICCg22cEAADoSsaY/dbajD85TjACAAD3CkpNTU0dUaezsNTbg5K1VkeOHNH27dsVERGh+fPna+jQoR2PlZaWKjs7W8OHD9fChQs7HgMAAOjtCEYAAOChtbS0/EnUuTMwdVdQ+u5r3rhxo8uC0tmzZ5WVlSVJWrRokUaNGtXxWFVVlbZu3aqmpiYtXrxYcXFxD/YHBwAA0MMRjAAAgMf0xqBUV1ennJwcXbx4UfPnz1dqamrHDDdv3tS2bdt06tQpzZs3TxMnTmRPEQAA6JMIRgAAwDE9LSjV1NTIWqvk5GQtXbq0Y6F1S0uLCgsLtXfvXk2ePFmzZ8+Wv7+/p/5YAAAAHEcwAgAAPVZ3B6Vdu3YpJydHkjR+/Hilpqbq9u3b2rFjhyIjI7VgwQINHjzYU98uAABAj3G3YOTjxDAAAAB38vX1VXh4uMLDwzt9vLm5+U+WcpeVlT10UJo2bZpycnI0d+5cNTU16Te/+Y1aW1uVmJioGTNmKDQ0tJu+cwAAgJ6JYAQAAHo8Pz+/Lg1K391mlpubq4CAAC1evFgxMTE6evSoNm/eLOnbK4/S0tI0ZMiQbvs+AQAAegpuSQMAAH3e94NSTU2NTp06pREjRuj69eu6fv26kpOTlZaWpqioKFVVVam4uFjHjh1TaGhox21rQUFBTn8rAAAAXYodRgAAAHdRU1OjI0eO6MiRI2pra1NqamrH1UUVFRU6cuSIysrKFB0drfHjx2vcuHHy8/NzemwAAIBHRjACAAC4B2utLl26pOLiYh09elTBwcEdVxf5+/urrKxMxcXFqqys1NixYzV+/HiNHj1aXl5eTo8OAADwUAhGAAAAD8DlcunMmTMqLi5WWVmZIiIilJaWpnHjxqm1tVXHjh1TcXGxamtrlZKSorS0NEVERDzQs7UBAAA4jWAEAADwkFpaWlRWVqYjR47o7NmzSkhI6Li6qLa2tuN2Noll2QAAoHchGAEAAHSBhoYGHTt2TEeOHNG1a9dYlg0AAHo1ghEAAEAXY1k2AADo7QhGAAAAHsKybAAA0FsRjAAAALoBy7KB7nPu3DmdPn1agYGBGjBggAYMGOD2uq+vr9MjAkCPRzACAADoZizLBjzjzJkzys3N1ZkzZ+75vsHBwQoKCuo0KHX2NpEJQH9DMAIAAHAQy7KBR2Ot1enTp5Wbm6tz585JkoYNG6bJkyertbVVDQ0Namho0O3btztev3HjhlpbWx/4axGZAPQnBCMAAIAegmXZwP2z1urUqVPKzc3V+fPnJUlhYWF6/PHHNW7cuHvezulyudTY2NgRkToLS0QmAP0ZwQgAAKCHYVk2cHfWWpWXlysvL09VVVWSpKFDh2revHlKTk726N4vIhOA/oRgBAAA0IOxLBv4lrVWx48fV15eni5duiRJGjJkiObOnauUlJQeG0wfJDLdunVLN27cUFtb2wN/HSITgK5GMAIAAOglWJaN/sjlcqm0tFR5eXm6cuWKJCk0NFSZmZlKS0vrsaHoURCZAPQEBCMAAIBeiGXZ6OtcLpeOHj2q/Px8Xbt2TdZaDRo0SHPmzNGECRPk7e3t9Ig9yv1Epu/e9lRk+n5gIjIBvRvBCAAAoJdjWTb6EpfLpeLiYuXn5+vmzZtyuVwKCgrS7NmzNWnSJEJRF7pbZOrsaiYiE9D/EIwAAAD6CJZlozdra2vT4cOHVVBQoMbGRllr5ePjo1mzZik9PV0+Pj5OjwgRmYD+hGAEAADQB91rWfbRo0d15MgRlmXDca2trTp06JAKCgrkcrn03b9DZsyYoYyMDEJAH0BkAnonghEAAEAfx7Js9EQtLS06cOCACgsL5e3tLWutWlpa9Nhjj2nKlCncNtnPEZkA5xGMAAAA+hGWZcNpLS0tKioqUmFhoQIDAyVJt27d0vTp0zV16lT5+/s7PCF6Kyci072eWY7IhN6MYAQAANBPsSwb3am5uVnffPONdu/erZCQEElSbW2tpk2bpmnTpikgIMDhCdEfOR2ZOgtMRCb0FAQjAACAfo5l2fCkpqYm7du3T3v27NHQoUMlSdXV1crIyNBjjz3WcZUR0Fvcb2T67oXIhN6KYAQAAIAOLMtGV7l9+7b27t2rb775RsOHD5ckXbx4Uenp6ZoxY4YGDBjg8IRA9yEyoTciGAEAAKBTLMvGw2hoaNCePXtUVFSkyMhIGWNUWVmpSZMmaebMmezEAu7TD0WmzkLTw0amgQMHauDAgUQm/AmCEQAAAO6JZdm4l1u3bmn37t06cOCARo0aJWOMTp8+rbS0NM2aNUvBwcFOjwj0eUQmdCWCEQAAAB4Iy7Jxp/r6ehUWFurgwYMdu61OnjyplJQUzZ49u2PBNYCeqbsj0/08sxyRqWcgGAEAAOChsCy7f6urq9OuXbtUXFysxMREeXl56fjx40pKStLs2bMVGhrq9IgAPMTlcun27dud7l/ydGS6W2AiMnU9ghEAAAAeGcuy+48bN26ooKBAR48eVVJSkry9vXXs2DElJiZqzpw5Gjx4sNMjAuiBiEy9D8EIAAAAXYpl2X1TTU2NCgoKVFpaqtTUVHl7e6u4uFhjxozRnDlzNHToUKdHBNDHEJmcRTACAACAx7Asu/e7fv268vPzVVZWpgkTJsjb21uHDh1SbGysMjMzFR4efteP/fLgBf23rWWqqr2tiNBA/R+LE/XcpMhunB5Af+N0ZEpLS+szAZ1gBAAAgG7Bsuzepbq6Wvn5+Tpx4oQmT54sHx8f7d+/X9HR0crMzNTw4cN/8OO/PHhBf7n5iG63/PEfYoG+3vp/nh9PNALQo3RlZFq8eLGmT5/ezd+BZxCMAAAA0K1Ylt2zXblyRfn5+aqoqFBGRoZ8fHz0zTffKCIiQnPnztWIESPu6/PM/PvtulB7+0+OR4YGatdfPN7VYwNAt7pbZEpISFBwcLDT43UJghEAAAAcw7LsnuPSpUvKy8vTuXPnNHXqVPn4+Gjv3r0aPny45s6dq4iIiAf6fHF/8Tt19i8KI+n03z/dJTMDADznbsHIx4lhAAAA0L94eXkpPj5e8fHxbsuy/+M//qNjWfbKlSs7lmVv3rxZEsuyu1JVVZXy8vJ04cIFTZ8+XbGxsdq9e7eGDh2qZcuWKSoq6qE+b0RoYKdXGEWEBj7qyAAAB3GFEQAAABzDsmzPO3/+vPLy8nTp0iXNmDFDvr6+KiwsVEhIiObOnatRo0Y90udnhxEA9G7ckgYAAIAejWXZXevcuXPKy8tTdXW1ZsyYIT8/PxUUFCgoKEhz585VXFxcl30tniUNAHovghEAAAB6hc6WZaelpSk1NVV+fn4sy/4B1lqdPXtWubm5qq2t1axZs+Tn56f8/Hz5+flp3rx5io+PZzcUAKADwQgAAAC9zveXZUdGRnZcXcSy7D+y1ur06dPKzc1VfX29Zs2aJX9/f+Xl5cnb21tz587VmDFj+t2fCwDg3ro9GBlj/pukZyQ1SzolaaW1ttYYEyupVFJZ+7vusda+c6/PRzACAADo3+5cln327NmOZdmjR4/uWJZ95MgRSf1nWba1VidPnlReXp4aGxs1e/Zs+fv7Kzc3Vy6XS/PmzdPYsWMJRQCAu3IiGC2StN1a22qM+QdJstb+eXsw2mKtTX2Qz0cwAgAAwHf6+7Jsa63Ky8uVl5en1tZWt1DU3NysuXPnKikpiVAEALgnR29JM8YslfSitfY1ghEAAAC6Un9alm2t1fHjx5WXlydrrTIzMztCUUNDgzIzM5WSkkIoAgDcN6eD0W8l/cpa+2F7MDomqVxSnaT/01qbf6/PQTACAADAD+nLy7JdLpdKSkqUn58vb29vZWZmKiAgQDt37lRdXZ0yMzOVmpraK74XAEDP4pFgZIzJkTSik4f+2lr7Vfv7/LWkDEnPW2utMcZf0kBr7TVjTLqkLyWlWGvrOvn8ayWtlaSYmJj0s2fPPvSsAAAA6D/6yrJsl8ulo0ePKj8/X/7+/m6h6Pr168rMzFRaWhqhCADw0By5wsgY86akdyTNt9Y23OV9dkr63621P3j5EFcYAQAA4GH0xmXZbW1tOnLkiPLz8zVw4EDNmTOnIxRdvXpVs2fP1sSJE+Xt7e3onACA3s+JpddPSPonSZnW2qt3HA+XdN1a22aMiZeUL2m8tfb6D30+ghEAAAAeVU9flt3W1qbDhw8rPz9foaGhblcUVVVVafbs2Zo0aZJ8fHy6bSYAQN/mRDA6Kclf0rX2Q3uste8YY16Q9HeSWiW1Sfoba+1v7/X5CEYAAADoSj1pWXZra6sOHjyoXbt2KSwszO2KosrKSs2aNUvp6emEIgBAl3N06XVXIBgBAADAE5xclt3S0qIDBw5o165dGjFihObMmdPxrGdnzpzRzJkzlZGRIV9f3y74TgEA+FMEIwAAAOAeumtZdnNzs/bv36/CwkJFRka6haJTp07pscce09SpUz12RRMAAN8hGAEAAAAPwBPLspubm/XNN99o9+7diomJ6QhFeXl5Ki8v17Rp0zRt2jT5+/t3x7cIAADBCAAAAHhYj7osu6mpSfv27dOePXsUFxenOXPmyM/PT/n5+SotLdXUqVM1ffp0BQQEOPhdAgD6I4IRAAAA0AUeZFl2YmKibty4oQMHDmjMmDGaNWuW/P39lZ+fr2PHjikjI0OPPfaYAgMDnf62AAD9FMEIAAAA6EL3syz7q6++krVWixYtUnJysgoKCnTkyBFNmjRJM2fO1IABA5z+NgAA/RzBCAAAAPCQzpZlp6amKisrS62trTLGyOVyKT09XbNmzdLAgQOdHhkAAEl3D0Y+TgwDAAAA9CVeXl6Kj49XfHx8x7LsQ4cO6fbt2/L19VVKSoquXLmi69evy8eHX8EBAD2fl9MDAAAAAH2Jr6+vUlNT9frrr2vVqlUaMGCAgoKC9NZbbyk0NFQbNmxQdXW102MCAPCDCEYAAACAh0RFRWnNmjU6f/68fvWrX2nevHmaMWOG3n//fZWXlzs9HgAAd0UwAgAAADwoKChIK1asUFhYmNavX6+IiAi98sor2rJli/Ly8tRbdooCAPoXghEAAADgYV5eXlq8eLHmzZunf//3f1dtba1Wr16t8vJyff7552pubnZ6RAAA3BCMAAAAgG4yfvx4vfHGG9qxY4cKCwv1xhtvyM/PTxs3blRNTY3T4wEA0IFgBAAAAHSj4cOHa82aNbp27Zo+/vhjLViwQJMnT9bGjRtVUVHh9HgAAEgiGAEAAADdLjAwUMuXL1dMTIzee+89RUVF6YUXXtDmzZu1e/du9hoBABxHMAIAAAAc4OXlpccff1xPPvmkPv74Y9XU1Gj16tUqLi7Wl19+qZaWFqdHBAD0YwQjAAAAwEHjxo3TypUrtXv3bhUUFOiNN95QW1ubPvjgA9XV1Tk9HgCgnyIYAQAAAA4LCwvT6tWr1dDQoI8//liLFi1SUlKSNmzYoHPnzjk9HgCgHyIYAQAAAD2Av7+/li1bprFjx2rDhg2Kjo7WkiVL9Nlnn6moqMjp8QAA/QzBCAAAAOghjDGaPXu2lixZol//+te6du2aVq5cqb1792rLli1qa2tzekQAQD9BMAIAAAB6mDFjxmjVqlU6ePCg8vLy9Oabb6q+vl6bNm1SfX290+MBAPoBghEAAADQAw0ePFirVq2StVYfffSRFi9erLi4OK1fv15VVVVOjwcA6OMIRgAAAEAP5evrq6VLl2rChAnauHGjoqOj9cQTT+ijjz7S4cOHnR4PANCHEYwAAACAHswYo+nTp+vFF1/Ul19+qWvXrunNN99Ubm6utm7dKpfL5fSIAIA+iGAEAAAA9AKxsbFas2aNjh8/rtzcXL355pu6evWqPvzwQzU0NDg9HgCgjyEYAQAAAL1ESEiI3nrrLfn7+3fsNRoxYoTWr1+vy5cvOz0eAKAPIRgBAAAAvYiPj4+WLFmiadOm6YMPPtCoUaM0b948bdq0SSUlJU6PBwDoIwhGAAAAQC+Unp6u5cuX6/e//72uXbum1157TVlZWdq2bZustU6PBwDo5QhGAAAAQC8VFRWlNWvW6MyZM9q5c6dWrFihyspKffLJJ2psbHR6PABAL0YwAgAAAHqxgQMH6o033tDgwYP18ccfa/HixQoNDdWGDRtUXV3t9HgAgF6KYAQAAAD0ct7e3nryySeVmZmpDz/8UKNGjdKMGTP0/vvvq7y83OnxAAC9EMEIAAAA6CPS0tL0+uuvKycnR9XV1XrppZe0ZcsW5eXlsdcIAPBACEYAAABAHzJy5EitWbNGly9fVm5url599VWVl5fr888/V3Nzs9PjAQB6CYIRAAAA0McMGDBAr732miIiIvTpp59q8eLF8vPz08aNG1VTU+P0eACAXoBgBAAAAPRBXl5eWrBggRYtWqRPP/1Uo0aN0uTJk7Vx40ZVVFQ4PR4AoIcjGAEAAAB9WHJyst566y0VFBTo2rVrWrp0qTZv3qw9e/aw1wgAcFcEIwAAAKCPCw8P1+rVq1VXV6fc3Fy98sorOnz4sL766iu1tLQ4PR4AoAciGAEAAAD9QEBAgF5++WWNHj1an332mRYuXKjW1lZ98MEHqqurc3o8AEAPQzACAAAA+gljjDIzM/XMM89o8+bNGjVqlJKSkrRhwwadO3fO6fEAAD0IwQgAAADoZxISEvT222+rqKhI165d01NPPaXPPvtMRUVFTo8GAOghCEYAAABAPzRkyBCtWrVKLS0tysvL0/PPP6+9e/dqy5Ytamtrc3o8AIDDCEYAAABAP+Xn56cXXnhBqamp+uKLLzR//nzV19dr06ZNqq+vd3o8AICDCEYAAABAP2aM0YwZM7R06VJt2bJF0dHRio2N1fr161VVVeX0eAAAhxCMAAAAACg+Pl6rV6/WsWPHdP36dc2fP18fffSRDh8+7PRoAAAHEIwAAAAASJJCQ0O1cuVK+fj4aNeuXXrmmWeUm5urrVu3yuVyOT0eAKAbEYwAAAAAdPD19dWSJUuUkZGhLVu2aO7cubp69ao+/PBDNTQ0OD0eAKCbEIwAAAAAuDHGaMqUKXrppZeUk5OjqKgojRgxQuvXr9fly5edHg8A0A08FoyMMf/ZGHPBGHOo/eWpOx77S2PMSWNMmTFmsadmAAAAAPDwYmJitGbNGlVUVOj69euaMWOGNm3apJKSEqdHAwB4mKevMPpna+3E9pffS5IxJlnSK5JSJD0h6WfGGG8PzwEAAADgIQQHB+vNN99UcHCw9u7dq8WLFysrK0vbtm2Ttdbp8QAAHuLELWnPSvrUWttkrT0t6aSkqQ7MAQAAAOA+eHt76+mnn9asWbO0detWzZo1S5WVlfrkk0/U2Njo9HgAAA/wdDBaZ4wpNsb8whgzuP1YpKTKO97nfPsxAAAAAD3YxIkT9dprr6mgoEAREREKCQnRhg0bVF1d7fRoAIAu9kjByBiTY4w52snLs5J+Lmm0pImSLkr6x+8+rJNP1em1rMaYtcaYImNM0dWrVx9lVAAAAABdICIiQmvWrNGlS5dUW1urSZMm6f3331d5ebnTowEAutAjBSNr7QJrbWonL19Zay9ba9ustS5J6/XH287OS4q+49NESaq6y+d/z1qbYa3NCA8Pf5RRAQAAAHSRoKAgvf766xo2bJiKioo0d+5cbdmyRXl5eew1AoA+wpPPkjbyjjeXSjra/vrXkl4xxvgbY+IkJUja56k5AAAAAHQ9Ly8vLVq0SPPnz9fOnTs1bdo0lZeX6/PPP1dzc7PT4wEAHpEndxj9V2PMEWNMsaR5kv43SbLWHpP0maQSSX+Q9K61ts2DcwAAAADwkNTUVL3xxhvav3+/RowYIR8fH23cuFE1NTVOjwYAeASmt1wympGRYYuKipweAwAAAEAnbt++rS+++ELNzc2KiYnRgQMH9Pzzzys+Pt7p0QAAP8AYs99am/H9455+ljQAAAAA/UBgYKCWL1+uUaNG6fDhw5o+fbo2b96sPXv2sNcIAHohH6cHAAAAANA3GGM0b948jRw5Ur/97W+Vnp6uw4cP69KlS3r66afl6+vr9IgAgPvEFUYAAAAAutS4ceO0cuVKlZaWKjw8XE1NTfrggw9UV1fn9GgAgPtEMAIAAADQ5cLCwrR69Wq1trbq5s2bioyM1IYNG3Tu3DmnRwMA3AeCEQAAAACP8Pf317JlyzRu3DiVlpZqwoQJ+uyzz8ST2QBAz8cOIwAAAAAeY4zRrFmzNHLkSH3xxRdKTk7W3r17denSJT355JPy9vZ2ekQAQCe4wggAAACAx40ePVqrVq1SZWWlhg4dqpqaGm3atEn19fVOjwYA6ATBCAAAAEC3GDx4sN5++235+/vr1q1bCg0N1fr161VVVeX0aACA7yEYAQAAAOg2vr6+eu655zRp0iSdOnVKCQkJ+uijj3T48GGnRwMA3IFgBAAAAKBbGWM0bdo0LVu2TGVlZRo9erRyc3O1detWuVwup8cDAIhgBAAAAMAho0aN0po1a3T9+nUNGjRI58+f14cffqiGhganRwOAfo9gBAAAAMAxISEheuuttzRkyBDdvn1bfn5+Wr9+vS5fvuz0aADQrxGMAAAAADjKx8dHzzzzjB577DFVVlYqMjJSmzZtUklJidOjAUC/RTACAAAA0COkp6dr+fLlHdFo69at2rZtm6y1To8GAP0OwQgAAABAjxEVFaU1a9aoublZAwcOVHl5uT755BM1NjY6PRoA9CsEIwAAAAA9ysCBA7VixQpFR0erqalJzc3N2rBhg6qrq50eDQD6DYIRAAAAgB7H29tbTzzxhB5//HFdvXpVgwYN0vvvv6/y8nKnRwOAfoFgBAAAAKDHSktL04oVK3T9+nUNHTpUX3/9tfLz89lrBAAeRjACAAAA0KONGDFCa9eulZ+fnwICAnTw4EF9/vnnam5udno0AOizCEYAAAAAerzAwEC9+uqrSkpKUnNzs6qrq7Vx40bV1NQ4PRoA9EkEIwAAAAC9gpeXl+bPn6+nn35a9fX18vb21saNG1VRUeH0aADQ5xCMAAAAAPQqSUlJeuutt9Tc3KzAwEB9/vnn2rNnD3uNAKALEYwAAAAA9Drh4eFavXq1hg4dKh8fHxUUFOirr75SS0uL06MBQJ9AMAIAAADQKwUEBOjll19WRkaG2tradPr0aX3wwQeqq6tzejQA6PUIRgAAAAB6LWOM5syZo+eff16tra26deuW1q9fr3Pnzjk9GgD0agQjAAAAAL1eQkKCVq1aJX9/f1lr9dFHH6moqMjpsQCg1yIYAQAAAOgThgwZolWrVikuLk7GGGVlZWnLli1qa2tzejQA6HUIRgAAAAD6DD8/Pz3//PPKzMyUJBUXF2vTpk2qr693eDIA6F0IRgAAAAD6FGOMHnvsMS1fvlx+fn66ePGi1q9fr6qqKqdHA4Beg2AEAAAAoE+Ki4vTmjVrFB4ervr6er3//vs6fPiw02MBQK9AMAIAAADQZw0aNEgrV65UWlqaWltb9eWXX2rr1q1yuVxOjwYAPRrBCAAAAECf5uPjoyVLluipp56Sl5eX9uzZow8//FANDQ1OjwYAPRbBCAAAAECfZ4zRlClT9NZbb2ngwIE6ffq01q9fr8uXLzs9GgD0SAQjAAAAAP1GdHS01q5dq+joaNXW1uq9995TSUmJ02MBQI9DMAIAAADQrwQHB+vNN9/UlClT5HK59Otf/1rbt2+Xtdbp0QCgxyAYAQAAAOh3vL299dRTT+nZZ5+VJOXn5+vTTz9VY2Ojw5MBQM9AMAIAAADQb02cOFFr1qxRcHCwysvLtX79elVXVzs9FgA4jmAEAAAAoF+LiIjQO++8o7i4OF2/fl0/+9nPVF5e7vRYAOAoghEAAACAfm/AgAF6/fXXNWPGDFlr9cUXX8jlcjk9FgA4xsfpAQAAAACgJ/Dy8tLChQsVHR0tX19feXnx/9cB9F8EIwAAAAC4w7hx45weAQAcRzIHAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALjx2LOkGWN+JSmx/c1QSbXW2onGmFhJpZLK2h/bY619x1NzAAAAAAAA4MF4LBhZa1/+7nVjzD9KunHHw6estRM99bUBAAAAAADw8DwWjL5jjDGSXpL0uKe/FgAAAAAAAB5dd+wwmi3psrX2xB3H4owxB40xucaY2d0wAwAAAAAAAO7TI11hZIzJkTSik4f+2lr7VfvryyV9csdjFyXFWGuvGWPSJX1pjEmx1tZ18vnXSlorSTExMY8yKgAAAAAAAO7TIwUja+2CH3rcGOMj6XlJ6Xd8TJOkpvbX9xtjTkkaK6mok8//nqT3JCkjI8M+yqwAAAAAAAC4P56+JW2BpOPW2vPfHTDGhBtjvNtfj5eUIKnCw3MAAAAAAADgPnl66fUrcr8dTZLmSPo7Y0yrpDZJ71hrr3t4DgAAAAAAANwnjwYja+1bnRz7jaTfePLrAgAAAAAA4OF1x7OkAQAAAAAAoBchGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3BCMAAAAAAAA4IZgBAAAAAAAADcEIwAAAAAAALghGAEAAAAAAMANwQgAAAAAAABuCEYAAAAAAABwQzACAAAAAACAG4IRAAAAAAAA3DxSMDLGLDPGHDPGuIwxGd977C+NMSeNMWXGmMV3HE83xhxpf+x/GmPMo8wAAAAAAACArvWoVxgdlfS8pLw7DxpjkiW9IilF0hOSfmaM8W5/+OeS1kpKaH954hFnAAAAAAAAQBd6pGBkrS211pZ18tCzkj611jZZa09LOilpqjFmpKQQa+1ua62VtEnSc48yAwAAAAAAALqWp3YYRUqqvOPt8+3HIttf//5xAAAAAAAA9BA+93oHY0yOpBGdPPTX1tqv7vZhnRyzP3D8bl97rb69fU0xMTH3mBQAAAAAAABd4Z7ByFq74CE+73lJ0Xe8HSWpqv14VCfH7/a135P0niRlZGTcNSwBAAAAAACg63jqlrSvJb1ijPE3xsTp2+XW+6y1FyXdNMZMb392tDck3e0qJQAAAAAAADjgkYKRMWapMea8pMck/c4Ys1WSrLXHJH0mqUTSHyS9a61ta/+wn0jaoG8XYZ+S9B+PMgMAAAAAAAC6lvn2ycp6voyMDFtUVOT0GAAAAAAAAH2GMWa/tTbj+8c9dUsaAAAAAAAAeimCEQAAAAAAANwQjAAAAAAAAOCm1+wwMsZclXTW6Tn6qTBJ1U4PAfQTnG9A9+KcA7oP5xvQvTjncL9GWWvDv3+w1wQjOMcYU9TZAiwAXY/zDehenHNA9+F8A7oX5xweFbekAQAAAAAAwA3BCAAAAAAAAG4IRrgf7zk9ANCPcL4B3YtzDug+nG9A9+KcwyNhhxEAAAAAAADccIURAAAAAAAA3BCM0MEYs8wYc8wY4zLGZHzvsb80xpw0xpQZYxbfcTzdGHOk/bH/aYwx3T850DcYY55oP8dOGmP+wul5gN7OGPMLY8wVY8zRO44NMcZkG2NOtP938B2PdfqzDsC9GWOijTE7jDGl7b9P/i/txznnAA8wxgQYY/YZYw63n3N/236ccw5dhmCEOx2V9LykvDsPGmOSJb0iKUXSE5J+Zozxbn/455LWSkpof3mi26YF+pD2c+pfJT0pKVnS8vZzD8DD+0B/+nPpLyRts9YmSNrW/va9ftYBuLdWSf/JWpskabqkd9vPK845wDOaJD1urZ0gaaKkJ4wx08U5hy5EMEIHa22ptbask4eelfSptbbJWnta0klJU40xIyWFWGt322+XYW2S9Fz3TQz0KVMlnbTWVlhrmyV9qm/PPQAPyVqbJ+n69w4/K+mX7a//Un/8udXpz7rumBPoC6y1F621B9pfvympVFKkOOcAj7Dfqm9/07f9xYpzDl2IYIT7ESmp8o63z7cfi2x//fvHATy4u51nALrWcGvtRenbf+BKGtZ+nHMQ6CLGmFhJkyTtFecc4DHGGG9jzCFJVyRlW2s559ClfJweAN3LGJMjaUQnD/21tfaru31YJ8fsDxwH8OA4nwBncQ4CXcAYM1DSbyT9r9bauh9Yb8k5Bzwia22bpInGmFBJXxhjUn/g3Tnn8MAIRv2MtXbBQ3zYeUnRd7wdJamq/XhUJ8cBPLi7nWcAutZlY8xIa+3F9lurr7Qf5xwEHpExxlffxqKPrLWb2w9zzgEeZq2tNcbs1Le7iTjn0GW4JQ3342tJrxhj/I0xcfp2ufW+9kscbxpjprc/O9obku52lRKAH/aNpARjTJwxxk/fLiX82uGZgL7oa0lvtr/+pv74c6vTn3UOzAf0Su2/C26UVGqt/ac7HuKcAzzAGBPefmWRjDGBkhZIOi7OOXQhrjBCB2PMUkn/r6RwSb8zxhyy1i621h4zxnwmqUTfPgPGu+2XP0rST/Tts9AESvqP9hcAD8ha22qMWSdpqyRvSb+w1h5zeCygVzPGfCJprqQwY8x5SX8j6e8lfWaMWSXpnKRlknSPn3UA7m2mpBWSjrTvVJGkvxLnHOApIyX9sv2ZzrwkfWat3WKM2S3OOXQR8+2TWwEAAAAAAADf4pY0AAAAAAAAuCEYAQAAAAAAwA3BCAAAAAAAAG4IRgAAAAAAAHBDMAIAAAAAAIAbghEAAAAAAADcEIwAAAAAAADghmAEAAAAAAAAN/8/LVWheoTJfr0AAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "fig, ax = plt.subplots(1, 1, figsize=(20,8))\n", - "ax.grid(False)\n", - "zips = []\n", - "cell.diam3d = cell.diam3d * 5\n", - "for x, y in cell.get_idx_polygons(projection=('x', 'y')):\n", - " zips.append(list(zip(x, y)))\n", - "polycol = PolyCollection(zips,\n", - " edgecolors='k',\n", - " facecolors='w')\n", - "# ax.add_collection(polycol)\n", - " \n", - "zips = []\n", - "for x, y in cell.get_pt3d_polygons(projection=('x', 'y')):\n", - " zips.append(list(zip(x, y)))\n", - "polycol = PolyCollection(zips,\n", - " edgecolors='gray',\n", - " facecolors='none')\n", - "ax.add_collection(polycol)\n", - "\n", - "\n", - "ax.plot(cell.x[stim.idx].mean(), cell.y[stim.idx].mean(), 'o')\n", - "\n", - "ax.axis('equal')" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -321,7 +173,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -349,7 +201,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -379,7 +231,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -404,7 +256,7 @@ " colors = [plt.get_cmap(cmap)(norm(v)) for v in V_m]\n", " zips = []\n", " for i in range(V_m.size):\n", - " inds = cell_geometry._c_ind == i\n", + " inds = cell_geometry._compartment_index == i\n", " zips.append(create_polygon(cell_geometry.x[inds, ].flatten(),\n", " cell_geometry.y[inds, ].flatten(),\n", " cell_geometry.d[inds, ].flatten()))\n", @@ -417,7 +269,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -445,7 +297,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -454,7 +306,7 @@ "Text(0, 0.5, '$y$ ($\\\\mu$m)')" ] }, - "execution_count": 15, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, @@ -527,7 +379,7 @@ "Hence the dendritic section receiving input about 1/6 of the way from its root \n", "is treated as 3 separate line sources with homogeneous current density per length unit each, even if \n", "the diameter and direction varies with position due to the fact \n", - "that each compartment may be composed of multiple pt3d segments. \n", + "that each compartment may be defined by multiple pt3d segments. \n", "\n", "The parameter `nseg=3` to `custom_fun` above can be changed via `custom_fun_args=[{nseg: }]` to affect the number of compartments per section. \n", "\n", diff --git a/examples/Example_Neuron_swc.ipynb b/examples/Example_Neuron_swc.ipynb index 7e920ea..6b8fdaa 100644 --- a/examples/Example_Neuron_swc.ipynb +++ b/examples/Example_Neuron_swc.ipynb @@ -14,7 +14,10 @@ "Its morphology is defined in the file `single_cell.swc`. \n", "\n", "The example produces comparable output as the `Example_Arbor_swc.ipynb` and `Example_LFPy_swc.ipynb` notebooks. \n", - "For this a couple of modified `CellGeometry` and `LineSourcePotential` classes have been defined below, similar to `Example_LFPy_pt3d.ipynb` and `Example_LFPy_swc.ipynb`" + "\n", + "This example showcase the use of `lfpykit.special.CellGeometryNeuron` and `lfpykit.special.LineSourcePotential`, \n", + "two modified classes inherited from `lfpykit.CellGeometry` and `lfpykit.LineSourcePotential`, respectively. \n", + "The modifications allow using the full geometry detail from the `.swc` morphology file for computing the extracellular potential." ] }, { @@ -208,177 +211,43 @@ { "cell_type": "code", "execution_count": 8, - "id": "ae9609b2-3752-4b28-9b2e-6a8537ee0556", + "id": "cc38538f-60ae-42cc-b54f-288ef0c84df7", "metadata": {}, "outputs": [], "source": [ - "class Pt3dCellGeometry(lfpykit.CellGeometry):\n", - " '''\n", - " Class inherited from ``lfpykit.CellGeometry`` for easier forward-model\n", - " predictions from NEURON that keeps track of pt3d information\n", - " for each compartment.\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: ``Cell`` object\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.CellGeometry\n", - " '''\n", - " def __init__(self, cell):\n", - " self.x3d, self.y3d, self.z3d, self.diam3d, self.arc3d = \\\n", - " self._collect_pt3d(cell)\n", - "\n", - " x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4\n", - " c_ind = np.array([], dtype=int) # tracks which compartment owns segment\n", - " idx = 0\n", - " for i, sec in enumerate(cell.all):\n", - " xp = (self.arc3d[i] / self.arc3d[i][-1])\n", - " seg_l = 1 / sec.nseg / 2 \n", - " for seg in sec:\n", - " seg_x = np.array([seg.x - seg_l, seg.x + seg_l])\n", - " ind = (xp > seg_x.min()) & (xp < seg_x.max())\n", - " seg_x = np.sort(np.r_[seg_x, xp[ind]])\n", - " \n", - " # interpolate to get values for start-, pt3d-, and end-point for search segment:\n", - " x_tmp = np.interp(seg_x, xp, self.x3d[i])\n", - " y_tmp = np.interp(seg_x, xp, self.y3d[i])\n", - " z_tmp = np.interp(seg_x, xp, self.z3d[i])\n", - " d_tmp = np.interp(seg_x, xp, self.diam3d[i])\n", - " \n", - " # store\n", - " x = np.row_stack([x, np.c_[x_tmp[:-1], x_tmp[1:]]])\n", - " y = np.row_stack([y, np.c_[y_tmp[:-1], y_tmp[1:]]])\n", - " z = np.row_stack([z, np.c_[z_tmp[:-1], z_tmp[1:]]])\n", - " d = np.row_stack([d, np.c_[d_tmp[:-1], d_tmp[1:]]])\n", - " \n", - " c_ind = np.r_[c_ind, np.ones(seg_x.size - 1) * idx]\n", - " idx += 1 \n", - " \n", - " super().__init__(x=x, y=y, z=z, d=d)\n", - " self._c_ind = c_ind\n", - "\n", - " def _collect_pt3d(self, cell):\n", - " \"\"\"collect the pt3d info, for each section\n", - "\n", - " Returns\n", - " -------\n", - " x3d: list of ndarray\n", - " x-coordinates from h.x3d()\n", - " y3d: list of ndarray\n", - " y-coordinates from h.y3d()\n", - " z3d: list of ndarray\n", - " z-coordinates from h.z3d()\n", - " diam3d: list of ndarray\n", - " diameter from h.diam3d()\n", - " arc3d: list of ndarray\n", - " arclength from h.arc3d()\n", - " \"\"\"\n", - " x3d, y3d, z3d, diam3d, arc3d = [], [], [], [], []\n", - "\n", - " for sec in cell.all:\n", - " n3d = int(neuron.h.n3d(sec=sec))\n", - " x3d_i, y3d_i, z3d_i, diam3d_i, arc3d_i = np.zeros((5, n3d))\n", - " for i in range(n3d):\n", - " x3d_i[i] = neuron.h.x3d(i, sec=sec)\n", - " y3d_i[i] = neuron.h.y3d(i, sec=sec)\n", - " z3d_i[i] = neuron.h.z3d(i, sec=sec)\n", - " diam3d_i[i] = neuron.h.diam3d(i, sec=sec)\n", - " arc3d_i[i] = neuron.h.arc3d(i, sec=sec)\n", - "\n", - " x3d.append(x3d_i)\n", - " y3d.append(y3d_i)\n", - " z3d.append(z3d_i)\n", - " diam3d.append(diam3d_i)\n", - " arc3d.append(arc3d_i)\n", - " \n", - " return x3d, y3d, z3d, diam3d, arc3d " + "# instantiate cell and run simulation\n", + "cell = CellSimulation()\n", + "cell.simulate()" ] }, { - "cell_type": "code", - "execution_count": 9, - "id": "056ed2f9-37c1-4c48-b9a7-fa80c0a737d1", + "cell_type": "markdown", + "id": "592fa242-ebfe-472b-81d6-c9a3badf7a6f", "metadata": {}, - "outputs": [], "source": [ - "class Pt3dLineSourcePotential(lfpykit.LineSourcePotential):\n", - " '''subclass of ``lfpykit.LineSourcePotential`` modified for\n", - " instances of ``Pt3dCellGeometry``.\n", - " Each compartment may consist of several segments, and this implementation\n", - " accounts for their contributions normalized by surface area, that is,\n", - " we assume constant transmembrane current density per area across each compartment\n", - " and constant current source density per unit length per segment\n", - " (inherent in the line-source approximation).\n", - "\n", - " Parameters\n", - " ----------\n", - " cell: object\n", - " ``Pt3dCellGeometry`` instance or similar.\n", - " x: ndarray of floats\n", - " x-position of measurement sites (µm)\n", - " y: ndarray of floats\n", - " y-position of measurement sites (µm)\n", - " z: ndarray of floats\n", - " z-position of measurement sites (µm)\n", - " sigma: float > 0\n", - " scalar extracellular conductivity (S/m)\n", - "\n", - " See also\n", - " --------\n", - " lfpykit.LineSourcePotential\n", - " '''\n", - "\n", - " def __init__(self, **kwargs):\n", - " super().__init__(**kwargs)\n", - " self._get_transformation_matrix = super().get_transformation_matrix\n", - "\n", - " def get_transformation_matrix(self):\n", - " '''Get linear response matrix\n", - "\n", - " Returns\n", - " -------\n", - " response_matrix: ndarray\n", - " shape (n_coords, n_CVs) ndarray\n", - " '''\n", - " M_tmp = self._get_transformation_matrix()\n", - " n_comps = np.unique(self.cell._c_ind).size\n", - " M = np.zeros((self.x.size, n_comps))\n", - " for i in range(n_comps):\n", - " inds = self.cell._c_ind == i\n", - " M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] /\n", - " self.cell.area[inds].sum())\n", - "\n", - " return M" + "## compute extracellular potentials\n", + "Combining the use of `lfpykit.special.CellGeometryNeuron` to represent the coordinates of `pt3d`-points for each compartment and \n", + "`lfpykit.special.LineSourcePotential` compute the extracellular potential" ] }, { "cell_type": "code", - "execution_count": 10, - "id": "cc38538f-60ae-42cc-b54f-288ef0c84df7", + "execution_count": 9, + "id": "307cdefe-0c31-4817-b735-da0ae67738fe", "metadata": {}, "outputs": [], "source": [ - "# instantiate cell and run simulation\n", - "cell = CellSimulation()\n", - "cell.simulate()" + "# create ``CellGeometryLFPyPt3d`` instance\n", + "cell_geometry = lfpykit.special.CellGeometryNeuron(cell)" ] }, { "cell_type": "code", - "execution_count": 11, - "id": "7b862c17-4897-4308-ac9c-7a70285f4fb4", + "execution_count": 10, + "id": "b1033964-349b-4b82-9d2a-f30cafadf6ab", "metadata": {}, "outputs": [], "source": [ - "###############################################################################\n", - "# compute extracellular potential using pt3d information\n", - "###############################################################################\n", - "\n", - "# create ``Pt3dCellGeometry`` instance\n", - "cell_geometry = Pt3dCellGeometry(cell)\n", - "\n", "# locations where extracellular potential is predicted \n", "dx = 1\n", "dz = 1\n", @@ -390,11 +259,11 @@ "\n", "\n", "# LineSourcePotential object, get mapping\n", - "lsp_pt3d = Pt3dLineSourcePotential(cell=cell_geometry, \n", + "lsp = lfpykit.special.LineSourcePotential(cell=cell_geometry, \n", " x=X.flatten(), \n", " y=Y.flatten(), \n", " z=Z.flatten())\n", - "M = lsp_pt3d.get_transformation_matrix()\n", + "M = lsp.get_transformation_matrix()\n", "\n", "# Extracellular potential in x,z-plane coordinates \n", "V_e = M @ cell.I_m" @@ -411,7 +280,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "id": "86be136a-36d1-4203-9c80-12d1070ecdf7", "metadata": {}, "outputs": [], @@ -440,7 +309,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "id": "3437196e-ef11-4c69-966f-2f32576484f1", "metadata": {}, "outputs": [], @@ -465,7 +334,7 @@ " colors = [plt.get_cmap(cmap)(norm(v)) for v in V_m]\n", " zips = []\n", " for i in range(V_m.size):\n", - " inds = cell_geometry._c_ind == i\n", + " inds = cell_geometry._compartment_index == i\n", " zips.append(create_polygon(cell_geometry.x[inds, ].flatten(),\n", " cell_geometry.y[inds, ].flatten(),\n", " cell_geometry.d[inds, ].flatten()))\n", @@ -478,7 +347,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "id": "399ebee1-1324-4e5d-baf5-4e88a80c2c4b", "metadata": {}, "outputs": [], @@ -507,7 +376,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "id": "10d97302-966c-480e-a85c-90a8af901c80", "metadata": {}, "outputs": [], @@ -540,7 +409,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 15, "id": "33b372e5-f2bc-48f7-957c-e44095a63606", "metadata": {}, "outputs": [], @@ -569,7 +438,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 16, "id": "efc724d8-44fa-4663-9191-663586c58837", "metadata": {}, "outputs": [ @@ -579,7 +448,7 @@ "Text(0, 0.5, '$y$ ($\\\\mu$m)')" ] }, - "execution_count": 17, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, @@ -658,7 +527,7 @@ "Hence the dendritic section receiving input about 1/6 of the way from its root \n", "is treated as 3 separate line sources with homogeneous current density per length unit each, even if \n", "the diameter and direction varies with position due to the fact \n", - "that each compartment may be composed of multiple pt3d segments. \n", + "that each compartment may be defined by multiple pt3d points. \n", "\n", "The parameter `nseg=3` above can be changed via to affect the number of compartments per section. " ] diff --git a/lfpykit/__init__.py b/lfpykit/__init__.py index 61080d9..c8ff4a8 100644 --- a/lfpykit/__init__.py +++ b/lfpykit/__init__.py @@ -61,13 +61,15 @@ * models * eegmegcalc * lfpcalc + * special """ -from .version import version as __version__ +from .version import version as __version__ # noqa: F401 -from .cellgeometry import CellGeometry +from .cellgeometry import CellGeometry # noqa: F401 from .models import LinearModel, CurrentDipoleMoment, PointSourcePotential, \ LineSourcePotential, RecExtElectrode, RecMEAElectrode, \ OneSphereVolumeConductor, LaminarCurrentSourceDensity, \ - VolumetricCurrentSourceDensity -from . import eegmegcalc + VolumetricCurrentSourceDensity # noqa: F401 +from . import eegmegcalc # noqa: F401 +from . import special # noqa: F401 diff --git a/lfpykit/special.py b/lfpykit/special.py new file mode 100644 index 0000000..b2c6074 --- /dev/null +++ b/lfpykit/special.py @@ -0,0 +1,350 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Copyright (C) 2022 Computational Neuroscience Group, NMBU. + +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, either version 3 of the License, or +(at your option) any later version. +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. +""" +import numpy as np +from .cellgeometry import CellGeometry +from . import models + + +class CellGeometryArbor(CellGeometry): + ''' + Class inherited from ``lfpykit.CellGeometry`` for easier forward-model + predictions in Arbor that keeps track of arbor.segment information + for each CV, stored in the array ``CellGeometryArbor._compartment_index`` + + This class is usable with the forward models implemented in + ``lfpykit.special.CurrentDipoleMoment``, + ``lfpykit.special.PointSourcePotential``, + and ``lfpykit.special.LineSourcePotential``. + + Parameters + ---------- + p: object + ``arbor.place_pwlin`` object with 3-d locations and cables in a + morphology (cf. ``arbor.place_pwlin``) + cables: list + ``list`` of corresponding ``arbor.cable`` objects where transmembrane + currents are recorded (cf. ``arbor.cable_probe_total_current_cell``) + + See also + -------- + lfpykit.CellGeometry + ''' + + def __init__(self, p, cables): + x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4 + c_ind = np.array([], dtype=int) # tracks which CV owns segment + for i, m in enumerate(cables): + segs = p.segments([m]) + for j, seg in enumerate(segs): + x = np.row_stack([x, [seg.prox.x, seg.dist.x]]) + y = np.row_stack([y, [seg.prox.y, seg.dist.y]]) + z = np.row_stack([z, [seg.prox.z, seg.dist.z]]) + d = np.row_stack( + [d, [seg.prox.radius * 2, seg.dist.radius * 2]]) + c_ind = np.r_[c_ind, i] + + super().__init__(x=x, y=y, z=z, d=d) + self._compartment_index = c_ind + + +class CellGeometryNeuron(CellGeometry): + ''' + Class inherited from ``lfpykit.CellGeometry`` for easier forward-model + predictions from NEURON that keeps track of pt3d information for each + compartment, stored in the array ``CellGeometryNeuron._compartment_index`` + + This class is usable with the forward models implemented in + ``lfpykit.special.CurrentDipoleMoment``, + ``lfpykit.special.PointSourcePotential``, + and ``lfpykit.special.LineSourcePotential``. + + Parameters + ---------- + cell: object + ``Cell`` object instantiated using the ``neuron.h.Import3d_GUI`` + functionality + + See also + -------- + lfpykit.CellGeometry + ''' + + def __init__(self, cell): + self.x3d, self.y3d, self.z3d, self.diam3d, self.arc3d = \ + self._collect_pt3d(cell) + + x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4 + c_ind = np.array([], dtype=int) # tracks which comp. owns segment + idx = 0 + for i, sec in enumerate(cell.all): + xp = (self.arc3d[i] / self.arc3d[i][-1]) + seg_l = 1 / sec.nseg / 2 + for seg in sec: + seg_x = np.array([seg.x - seg_l, seg.x + seg_l]) + ind = (xp > seg_x.min()) & (xp < seg_x.max()) + seg_x = np.sort(np.r_[seg_x, xp[ind]]) + + # interpolate to get values for start-, pt3d-, and + # end-point for search segment: + x_tmp = np.interp(seg_x, xp, self.x3d[i]) + y_tmp = np.interp(seg_x, xp, self.y3d[i]) + z_tmp = np.interp(seg_x, xp, self.z3d[i]) + d_tmp = np.interp(seg_x, xp, self.diam3d[i]) + + # store + x = np.row_stack([x, np.c_[x_tmp[:-1], x_tmp[1:]]]) + y = np.row_stack([y, np.c_[y_tmp[:-1], y_tmp[1:]]]) + z = np.row_stack([z, np.c_[z_tmp[:-1], z_tmp[1:]]]) + d = np.row_stack([d, np.c_[d_tmp[:-1], d_tmp[1:]]]) + + c_ind = np.r_[c_ind, np.ones(seg_x.size - 1) * idx] + idx += 1 + + super().__init__(x=x, y=y, z=z, d=d) + self._compartment_index = c_ind + + def _collect_pt3d(self, cell): + """collect the pt3d info, for each section + + Returns + ------- + x3d: list of ndarray + x-coordinates from h.x3d() + y3d: list of ndarray + y-coordinates from h.y3d() + z3d: list of ndarray + z-coordinates from h.z3d() + diam3d: list of ndarray + diameter from h.diam3d() + arc3d: list of ndarray + arclength from h.arc3d() + """ + import neuron + + x3d, y3d, z3d, diam3d, arc3d = [], [], [], [], [] + + for sec in cell.all: + n3d = int(neuron.h.n3d(sec=sec)) + x3d_i, y3d_i, z3d_i, diam3d_i, arc3d_i = np.zeros((5, n3d)) + for i in range(n3d): + x3d_i[i] = neuron.h.x3d(i, sec=sec) + y3d_i[i] = neuron.h.y3d(i, sec=sec) + z3d_i[i] = neuron.h.z3d(i, sec=sec) + diam3d_i[i] = neuron.h.diam3d(i, sec=sec) + arc3d_i[i] = neuron.h.arc3d(i, sec=sec) + + x3d.append(x3d_i) + y3d.append(y3d_i) + z3d.append(z3d_i) + diam3d.append(diam3d_i) + arc3d.append(arc3d_i) + + return x3d, y3d, z3d, diam3d, arc3d + + +class CellGeometryLFPyPt3d(CellGeometry): + ''' + Class inherited from ``lfpykit.CellGeometry`` for easier forward-model + predictions in LFPy that keeps track of pt3d information + for each compartment. ``LFPy.Cell`` must be instantiated with the + keyword argument ``pt3d=True``. + + This class is usable with the forward models implemented in + ``lfpykit.special.CurrentDipoleMoment``, + ``lfpykit.special.PointSourcePotential``, + and ``lfpykit.special.LineSourcePotential``. + + Parameters + ---------- + cell: object + ``LFPy.Cell`` like object + + See also + -------- + lfpykit.CellGeometry + ''' + + def __init__(self, cell): + x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4 + c_ind = np.array([], dtype=int) # tracks which comp. owns segment + idx = 0 + for i, sec in enumerate(cell.allseclist): + xp = (cell.arc3d[i] / cell.arc3d[i][-1]) + seg_l = 1 / sec.nseg / 2 + for seg in sec: + seg_x = np.array([seg.x - seg_l, seg.x + seg_l]) + ind = (xp > seg_x.min()) & (xp < seg_x.max()) + seg_x = np.sort(np.r_[seg_x, xp[ind]]) + + # interpolate to get values for start-, pt3d-, and + # end-point for search segment: + x_tmp = np.interp(seg_x, xp, cell.x3d[i]) + y_tmp = np.interp(seg_x, xp, cell.y3d[i]) + z_tmp = np.interp(seg_x, xp, cell.z3d[i]) + d_tmp = np.interp(seg_x, xp, cell.diam3d[i]) + + # store + x = np.row_stack([x, np.c_[x_tmp[:-1], x_tmp[1:]]]) + y = np.row_stack([y, np.c_[y_tmp[:-1], y_tmp[1:]]]) + z = np.row_stack([z, np.c_[z_tmp[:-1], z_tmp[1:]]]) + d = np.row_stack([d, np.c_[d_tmp[:-1], d_tmp[1:]]]) + + c_ind = np.r_[c_ind, np.ones(seg_x.size - 1) * idx] + idx += 1 + + super().__init__(x=x, y=y, z=z, d=d) + self._compartment_index = c_ind + + +class CurrentDipoleMoment(models.CurrentDipoleMoment): + '''subclass of ``lfpykit.CurrentDipoleMoment`` modified for + instances of ``lfpykit.special.CellGeometry*``. + Each compartment may consist of several segments, and this implementation + accounts for their contributions normalized by surface area, that is, + we assume constant transmembrane current density per area across each + compartment and constant current source density per area per segment. + + Parameters + ---------- + cell: object + ``special.CellGeometry*`` instance or similar. + x: ndarray of floats + x-position of measurement sites (µm) + y: ndarray of floats + y-position of measurement sites (µm) + z: ndarray of floats + z-position of measurement sites (µm) + sigma: float > 0 + scalar extracellular conductivity (S/m) + + See also + -------- + lfpykit.CurrentDipoleMoment + ''' + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def get_transformation_matrix(self): + '''Get linear response matrix + + Returns + ------- + response_matrix: ndarray + shape (n_coords, n_compartments) ndarray + ''' + M_tmp = super().get_transformation_matrix() + n_compartments = np.unique(self.cell._compartment_index).size + M = np.zeros((self.x.size, n_compartments)) + for i in range(n_compartments): + inds = self.cell._compartment_index == i + M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] / + self.cell.area[inds].sum()) + return M + + +class LineSourcePotential(models.LineSourcePotential): + '''subclass of ``lfpykit.LineSourcePotential`` modified for + instances of ``lfpykit.special.CellGeometry*``. + Each compartment may consist of several segments, and this implementation + accounts for their contributions normalized by surface area, that is, + we assume constant transmembrane current density per area across each + compartment and constant current source density per unit length per segment + (inherent in the line-source approximation). + + Parameters + ---------- + cell: object + ``special.CellGeometry*`` instance or similar. + x: ndarray of floats + x-position of measurement sites (µm) + y: ndarray of floats + y-position of measurement sites (µm) + z: ndarray of floats + z-position of measurement sites (µm) + sigma: float > 0 + scalar extracellular conductivity (S/m) + + See also + -------- + lfpykit.LineSourcePotential + ''' + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def get_transformation_matrix(self): + '''Get linear response matrix + + Returns + ------- + response_matrix: ndarray + shape (n_coords, n_compartments) ndarray + ''' + M_tmp = super().get_transformation_matrix() + n_compartments = np.unique(self.cell._compartment_index).size + M = np.zeros((self.x.size, n_compartments)) + for i in range(n_compartments): + inds = self.cell._compartment_index == i + M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] / + self.cell.area[inds].sum()) + return M + + +class PointSourcePotential(models.PointSourcePotential): + '''subclass of ``lfpykit.PointSourcePotential`` modified for + instances of ``lfpykit.special.CellGeometry*``. + Each compartment may consist of several segments, and this implementation + accounts for their contributions normalized by surface area, that is, + we assume constant transmembrane current density per area across each + compartment and constant current source density per area per segment. + + Parameters + ---------- + cell: object + ``lfpykit.special.CellGeometry*`` instance or similar. + x: ndarray of floats + x-position of measurement sites (µm) + y: ndarray of floats + y-position of measurement sites (µm) + z: ndarray of floats + z-position of measurement sites (µm) + sigma: float > 0 + scalar extracellular conductivity (S/m) + + See also + -------- + lfpykit.PointSourcePotential + ''' + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def get_transformation_matrix(self): + '''Get linear response matrix + + Returns + ------- + response_matrix: ndarray + shape (n_coords, n_compartments) ndarray + ''' + M_tmp = super().get_transformation_matrix() + n_compartments = np.unique(self.cell._compartment_index).size + M = np.zeros((self.x.size, n_compartments)) + for i in range(n_compartments): + inds = self.cell._compartment_index == i + M[:, i] = M_tmp[:, inds] @ (self.cell.area[inds] / + self.cell.area[inds].sum()) + return M