diff --git a/CMakeLists.txt b/CMakeLists.txt index 4ab316c1c5..6fd242ad28 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,6 +45,23 @@ set(NMODL_ClangFormat_DEPENDENCIES pyastgen parser-gen FORCE) add_subdirectory(cmake/hpc-coding-conventions/cpp) +# ============================================================================= +# Format & execute ipynb notebooks in place (pip install nbconvert clean-ipynb) +# ============================================================================= +add_custom_target(nb-format + jupyter + nbconvert + --to + notebook + --execute + --inplace + --ExecutePreprocessor.timeout=360 + "${CMAKE_SOURCE_DIR}/docs/notebooks/*.ipynb" + && + clean_ipynb + --keep-output + "${CMAKE_SOURCE_DIR}/docs/notebooks/*.ipynb") + # ============================================================================= # Include cmake modules # ============================================================================= diff --git a/docs/index.rst b/docs/index.rst index 2db6af6b5d..063a9ef719 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -22,9 +22,13 @@ Welcome to nmodl's documentation! :maxdepth: 3 :caption: Jupyter Notebooks: - notebooks/nmodl-python-sympy-examples - notebooks/nmodl-python-tutorial - notebooks/kinetic-schemes + notebooks/nmodl-python-tutorial.ipynb + notebooks/nmodl-odes-overview.ipynb + notebooks/nmodl-kinetic-schemes.ipynb + notebooks/nmodl-sympy-solver.ipynb + notebooks/nmodl-linear-solver.ipynb + notebooks/nmodl-nonlinear-solver.ipynb + notebooks/nmodl-sympy-conductance.ipynb .. toctree:: :maxdepth: 2 diff --git a/docs/notebooks/README.md b/docs/notebooks/README.md new file mode 100644 index 0000000000..ffd4f708e8 --- /dev/null +++ b/docs/notebooks/README.md @@ -0,0 +1,14 @@ +## NMODL jupyter notebooks + +To get started with the NMODL python interface: + - [nmodl-python-tutorial.ipynb](nmodl-python-tutorial.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-python-tutorial.ipynb) + +For an overview of ODEs in NODL: + - [nmodl-odes-overview.ipynb](nmodl-odes-overview.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-odes-overview.ipynb) + +For more specific implementation details: + - [nmodl-kinetic-schemes.ipynb](nmodl-kinetic-schemes.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-kinetic-schemes.ipynb) + - [nmodl-sympy-solver.ipynb](nmodl-sympy-solver.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-sympy-solver.ipynb) + - [nmodl-linear-solver.ipynb](nmodl-linear-solver.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-linear-solver.ipynb) + - [nmodl-nonlinear-solver.ipynb](nmodl-nonlinear-solver.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-nonlinear-solver.ipynb) + - [nmodl-sympy-conductance.ipynb](nmodl-sympy-conductance.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/BlueBrain/nmodl/blob/add_sympy_documentation/docs/notebooks/nmodl-sympy-conductance.ipynb) diff --git a/docs/notebooks/kinetic-schemes.ipynb b/docs/notebooks/kinetic-schemes.ipynb deleted file mode 100644 index db2dbd6090..0000000000 --- a/docs/notebooks/kinetic-schemes.ipynb +++ /dev/null @@ -1,243 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### NMODL Kinetic Scheme\n", - "\n", - "Some definitions of reaction kinetics & mass action laws that are relevant for how `KINETIC` blocks are translated into a system of ODEs in NMODL" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Reaction Kinetics\n", - "\n", - "- We consider a set of reaction species $A_j$, with corresponding Molar concentrations $Y_j$\n", - "- They react according to a set of reaction equations:\n", - "\n", - "$$\\sum_j \\nu_{ij}^L A_j \\overset{k_i}{\\rightarrow} \\sum_j \\nu_{ij}^R A_j$$\n", - "\n", - "where\n", - "- $k_i$ is the rate coefficient\n", - "- $\\nu_{ij}^L$, $\\nu_{ij}^R$ are stochiometric coefficients - must be positive integers (including zero)\n", - "\n", - "***\n", - "#### Law of Mass Action\n", - "\n", - "- This allows us to convert these reaction equations to a set of ODEs\n", - "\n", - "$$\\frac{dY_j}{dt} = \\sum_i \\Delta \\nu_{ij} r_i$$\n", - "\n", - "where $\\Delta \\nu_{ij} = \\nu_{ij}^R - \\nu_{ij}^L$, and\n", - "$$r_i = k_i \\prod_j Y_j^{\\nu_{ij}^R}$$\n", - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### KINETIC block format\n", - "\n", - "A reaction equation is specifed in the mod file as\n", - "```\n", - "~ A0 + 3A1 + 2A2 + ... <-> 2A0 + A1 + ... (kf, kb)\n", - "```\n", - "where\n", - "- `A0` etc are the species $A_j$\n", - "- the integer preceeding a species (with or without a space) is the corresponding stochiometric coefficient $\\nu_{ij}$ (by default 1 if not present)\n", - "- `kf` is the forwards reaction rate $k^{(f)}_j$\n", - "- `kb` is the backwards reaction rate $k^{(b)}_j$, i.e. the reaction rate for the same reaction with the LHS and RHS exchanged\n", - "\n", - "***\n", - "We can convert these statements to a system of ODEs:\n", - "$$\\frac{dY_j}{dt} = \\sum_i \\Delta \\nu_{ij} (r^{(f)}_i - r^{(b)}_i)$$\n", - "\n", - "where $\\Delta \\nu_{ij} = \\nu_{ij}^R - \\nu_{ij}^L$, and\n", - "$$\n", - "\\begin{aligned}\n", - "r^{(f)}_i &= k^{(f)}_i \\prod_j Y_j^{\\nu_{ij}^{L}}\\\\\n", - "r^{(b)}_i &= k^{(b)}_i \\prod_j Y_j^{\\nu_{ij}^{R}}\n", - "\\end{aligned}\n", - "$$\n", - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Example\n", - "\n", - "Given the following KINETIC statement\n", - "```\n", - "~ h <-> m (a,b)\n", - "```\n", - "We have 2 state variables, and 1 reaction statement, so \n", - "- $i = \\{0\\}$\n", - "- $j = \\{0, 1\\}$\n", - "\n", - "i.e. \n", - "\n", - "- the state vector $Y$ has 2 elements\n", - "$$\n", - "Y = \\left(\n", - "\\begin{aligned}\n", - "h \\\\\n", - "m \n", - "\\end{aligned}\n", - "\\right)\n", - "$$\n", - "\n", - "- the stochiometric coefficients are 1x2 matrices\n", - "$$\n", - "\\nu^L = \\left(\n", - "\\begin{aligned}\n", - "1 && 0\n", - "\\end{aligned}\n", - "\\right)\n", - "$$\n", - "$$\n", - "\\nu^R = \\left(\n", - "\\begin{aligned}\n", - "0 && 1\n", - "\\end{aligned}\n", - "\\right)\n", - "$$\n", - "\n", - "- the rate vectors contain 1 element:\n", - "$$\n", - "k^{(f)} = a\n", - "$$\n", - "$$\n", - "k^{(b)} = b\n", - "$$\n", - "\n", - "Using these we can construct the forwards and blackwards fluxes:\n", - "$$\n", - "r^{(f)} = a h\n", - "$$\n", - "$$\n", - "r^{(b)} = b m\n", - "$$\n", - "\n", - "and finally we find the ODEs in matrix form:\n", - "$$\n", - "\\frac{dY}{dt} =\n", - "\\left(\n", - "\\begin{aligned}\n", - "-1 && 1\n", - "\\end{aligned}\n", - "\\right)\n", - "(ah - bm)\n", - "$$\n", - "\n", - "which in terms of the state variables can be written:\n", - "$$\n", - "\\begin{aligned}\n", - "\\frac{dh}{dt} &= bm - ah \\\\\n", - "\\frac{dm}{dt} &= ah - bm\n", - "\\end{aligned}\n", - "$$\n", - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Other types of reaction statement\n", - "\n", - "Can also have a reaction statement of the form\n", - "```\n", - "~ h << (a)\n", - "```\n", - "where the LHS must be a state variable, and the RHS is an expression inside parentheses.\n", - "\n", - "The meaning of this statement is to add `a` to the differential equation for `h`, i.e.\n", - "$$\n", - "\\frac{dh}{dt} += a\n", - "$$\n", - "\n", - "***\n", - "Finally there is a statement of the form\n", - "```\n", - "x + 2y + ... -> (a)\n", - "```\n", - "which is a one-way reaction statement with no reaction products.\n", - "\n", - "This is just syntactic sugar for a special case of the standard `<->` reaction equation, where the backwards rate is set to zero and there are no states on the RHS of the reaction, so the above is equivalent to\n", - "```\n", - "x + 2y + ... <-> (a, 0)\n", - "```\n", - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### [TODO] CONSERVE\n", - "\n", - "Converse statement allows specification of a conservation law, e.g.\n", - "\n", - "```\n", - "CONSERVE h + m = 0\n", - "```\n", - "This can then be used to eliminate one of these variables from the set of ODEs, for example by removing the dmdt equation, and instead calculating m in terms of the algebraic relation above.\n", - "\n", - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### [TODO] fflux / bflux\n", - "\n", - "Directly after each kinetic statement, NEURON provides the fflux and bflux variable that can be referenced inside the mod file (!). Note that it always refers to the fflux/bflux from the preceeding kinetic statement.\n", - "\n", - "So we should parse all non-kinetic statements in the mod file following a kinetic statement (until we get to a new kinetic statement), and replace any fflux or bflux variables with corresponding expression.\n", - "\n", - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### [TODO] inlining\n", - "\n", - "The usual issues about inlining function calls, e.g. when the rate is a function call.\n", - "\n", - "In principle also allowed for it to be a function of state variables, in which case vital to inline this information for the ODEs (although we could also just say that this is not supported, as not necessary for the user to write the mod file in this way, state variable dependence can always be written explicitly).\n", - "\n", - "***" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.7" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/notebooks/nmodl-kinetic-schemes.ipynb b/docs/notebooks/nmodl-kinetic-schemes.ipynb new file mode 100644 index 0000000000..a319147bea --- /dev/null +++ b/docs/notebooks/nmodl-kinetic-schemes.ipynb @@ -0,0 +1,449 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### NMODL Kinetic Scheme\n", + "\n", + "This notebook describes the reaction kinetics & mass action laws that apply within `KINETIC` blocks, and the implementation of the `KineticBlockVisitor` in NMODL which transforms `KINETIC` blocks into `DERIVATIVE` blocks containing an equivalent system of ODEs.\n", + "\n", + "For a higher level overview of the approach to solving ODEs in NMODL, please see the [nmodl-odes-overview](nmodl-odes-overview.ipynb) notebook. \n", + "\n", + "For a more general tutorial on using the NMODL python interface, please see the [tutorial notebook](nmodl-python-tutorial.ipynb).\n", + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Reaction Kinetics\n", + "\n", + "- We consider a set of reaction species $A_j$, with corresponding Molar concentrations $Y_j$\n", + "- They react according to a set of reaction equations:\n", + "\n", + "$$\\sum_j \\nu_{ij}^L A_j \\overset{k_i}{\\rightarrow} \\sum_j \\nu_{ij}^R A_j$$\n", + "\n", + "where\n", + "- $k_i$ is the rate coefficient\n", + "- $\\nu_{ij}^L$, $\\nu_{ij}^R$ are stoichiometric coefficients - must be positive integers (including zero)\n", + "***\n", + "#### Law of Mass Action\n", + "\n", + "- This allows us to convert these reaction equations to a set of ODEs\n", + "\n", + "$$\\frac{dY_j}{dt} = \\sum_i \\Delta \\nu_{ij} r_i$$\n", + "\n", + "where $\\Delta \\nu_{ij} = \\nu_{ij}^R - \\nu_{ij}^L$, and\n", + "$$r_i = k_i \\prod_j Y_j^{\\nu_{ij}^R}$$\n", + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### KINETIC block format\n", + "\n", + "A reaction equation is specifed in the mod file as\n", + "```\n", + "~ A0 + 3A1 + 2A2 + ... <-> 2A0 + A1 + ... (kf, kb)\n", + "```\n", + "where\n", + "- `A0` etc are the species $A_j$\n", + "- the integer preceeding a species (with or without a space) is the corresponding stochiometric coefficient $\\nu_{ij}$ (implicitly 1 if not specified)\n", + "- `kf` is the forwards reaction rate $k^{(f)}_j$\n", + "- `kb` is the backwards reaction rate $k^{(b)}_j$, i.e. the reaction rate for the same reaction with the LHS and RHS exchanged\n", + "***\n", + "We can convert these statements to a system of ODEs using the law of Mass Action:\n", + "$$\\frac{dY_j}{dt} = \\sum_i \\Delta \\nu_{ij} (r^{(f)}_i - r^{(b)}_i)$$\n", + "\n", + "where $\\Delta \\nu_{ij} = \\nu_{ij}^R - \\nu_{ij}^L$, and\n", + "$$\n", + "\\begin{align}\n", + "r^{(f)}_i &= k^{(f)}_i \\prod_j Y_j^{\\nu_{ij}^{L}} \\\\\n", + "r^{(b)}_i &= k^{(b)}_i \\prod_j Y_j^{\\nu_{ij}^{R}}\n", + "\\end{align}\n", + "$$\n", + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Other types of reaction statement\n", + "There is also have a reaction statement of the form\n", + "```\n", + "~ h << (a)\n", + "```\n", + "where the LHS must be a state variable, and the RHS is an expression inside parentheses.\n", + "\n", + "The meaning of this statement is to add `a` to the differential equation for `h`, i.e.\n", + "$$\n", + "\\frac{dh}{dt} += a\n", + "$$\n", + "***\n", + "Finally there is a statement of the form\n", + "```\n", + "~ x + 2y + ... -> (a)\n", + "```\n", + "which is a one-way reaction statement with no reaction products.\n", + "\n", + "This is just syntactic sugar for a special case of the standard `<->` reaction equation, where the backwards rate is set to zero and there are no states on the RHS of the reaction, so the above is equivalent to\n", + "```\n", + "~ x + 2y + ... <-> (a, 0)\n", + "```\n", + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### f_flux / b_flux variables\n", + "\n", + "Within the `KINETIC` block in the MOD file, the user can make use of the `f_flux` and `b_flux` variables, which refer to the forwards $r^{(f)}$ and backwards $r^{(b)}$ fluxes of the preceeding reaction statement.\n", + "\n", + "The `KineticBlockVisitor` substitutes the current expression for these fluxes for these variables within the `KINETIC` block.\n", + "\n", + "If these variables are referenced before a reaction statement then they are assumed to be zero.\n", + "***\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### [TODO] CONSERVE\n", + "Converse statement allows specification of a conservation law, e.g.\n", + "\n", + "```\n", + "CONSERVE h + m + z = 1\n", + "```\n", + "In NEURON, the ODE for the last state variable on the rhs of this expression is replaced with this algebraic expression, so in this case instead of replacing $z' = ...$ with the forwards or backwards Euler equation, it would be replaced with $z = 1 - h - m$\n", + "\n", + "In order to be consistent with NEURON, in particular the way `STEADYSTATE` is implemented, we should do this in the same way.\n", + "***\n", + "#### Implementation Tests\n", + "\n", + " - The unit tests may be helpful to understand what these functions are doing\n", + " - `KineticBlockVisitor` tests are located in [test/visitor/visitor.cpp](https://github.com/BlueBrain/nmodl/blob/master/test/visitor/visitor.cpp#L2128)\n", + " ***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Examples" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import nmodl.dsl as nmodl\n", + "\n", + "\n", + "def run_kinetic_visitor_and_return_derivative(mod_string):\n", + " # parse NMDOL file (supplied as a string) into AST\n", + " driver = nmodl.NmodlDriver()\n", + " driver.parse_string(mod_string)\n", + " AST = driver.ast()\n", + " # run SymtabVisitor to generate Symbol Table\n", + " nmodl.symtab.SymtabVisitor().visit_program(AST)\n", + " # constant folding, inlining & local variable renaming passes\n", + " nmodl.visitor.ConstantFolderVisitor().visit_program(AST)\n", + " nmodl.visitor.InlineVisitor().visit_program(AST)\n", + " nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)\n", + " # run KINETIC block visitor\n", + " nmodl.visitor.KineticBlockVisitor().visit_program(AST)\n", + " # return new DERIVATIVE block\n", + " return nmodl.to_nmodl(\n", + " nmodl.visitor.AstLookupVisitor().lookup(\n", + " AST, nmodl.ast.AstNodeType.DERIVATIVE_BLOCK\n", + " )[0]\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 1\n", + "Given the following KINETIC statement\n", + "```\n", + "~ h <-> m (a,b)\n", + "```\n", + "We have 2 state variables, and 1 reaction statement, so \n", + "- $i = \\{0\\}$\n", + "- $j = \\{0, 1\\}$\n", + "\n", + "i.e. \n", + "- the state vector $Y$ has 2 elements\n", + "$$\n", + "Y = \\left(\n", + "\\begin{align}\n", + "h \\\\\n", + "m \n", + "\\end{align}\n", + "\\right)\n", + "$$\n", + "- the stoichiometric coefficients are 1x2 matrices\n", + "$$\n", + "\\nu^L = \\left(\n", + "\\begin{align}\n", + "1 && 0\n", + "\\end{align}\n", + "\\right)\n", + "$$\n", + "$$\n", + "\\nu^R = \\left(\n", + "\\begin{align}\n", + "0 && 1\n", + "\\end{align}\n", + "\\right)\n", + "$$\n", + "- the rate vectors contain 1 element:\n", + "$$\n", + "k^{(f)} = a\n", + "$$\n", + "$$\n", + "k^{(b)} = b\n", + "$$\n", + "\n", + "Using these we can construct the forwards and blackwards fluxes:\n", + "$$\n", + "r^{(f)} = a h\n", + "$$\n", + "$$\n", + "r^{(b)} = b m\n", + "$$\n", + "and finally we find the ODEs in matrix form:\n", + "$$\n", + "\\frac{dY}{dt} =\n", + "\\left(\n", + "\\begin{align}\n", + "-1 && 1\n", + "\\end{align}\n", + "\\right)\n", + "(ah - bm)\n", + "$$\n", + "which in terms of the state variables can be written:\n", + "$$\n", + "\\begin{align}\n", + "\\frac{dh}{dt} &= bm - ah \\\\\n", + "\\frac{dm}{dt} &= ah - bm\n", + "\\end{align}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DERIVATIVE kin {\n", + " h' = (-1*(a*h-b*m))\n", + " m' = (1*(a*h-b*m))\n", + "}\n" + ] + } + ], + "source": [ + "ex1 = \"\"\"\n", + "STATE {\n", + " h m\n", + "}\n", + "KINETIC kin {\n", + " ~ h <-> m (a,b)\n", + "}\n", + "\"\"\"\n", + "print(run_kinetic_visitor_and_return_derivative(ex1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 2\n", + "Annihilation reaction statement\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DERIVATIVE kin {\n", + " x' = (-1*(a*x))\n", + "}\n" + ] + } + ], + "source": [ + "ex2 = \"\"\"\n", + "STATE {\n", + " x\n", + "}\n", + "KINETIC kin {\n", + " ~ x -> (a)\n", + "}\n", + "\"\"\"\n", + "print(run_kinetic_visitor_and_return_derivative(ex2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 3\n", + "`<<` reaction statement" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DERIVATIVE kin {\n", + " x' = (a)\n", + "}\n" + ] + } + ], + "source": [ + "ex3 = \"\"\"\n", + "STATE {\n", + " x\n", + "}\n", + "KINETIC kin {\n", + " ~ x << (a)\n", + "}\n", + "\"\"\"\n", + "print(run_kinetic_visitor_and_return_derivative(ex3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 4\n", + "Annihilation & `<<` reaction statement for the same state variable" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DERIVATIVE kin {\n", + " x' = (a)+(-1*(b*x))\n", + "}\n" + ] + } + ], + "source": [ + "ex4 = \"\"\"\n", + "STATE {\n", + " x\n", + "}\n", + "KINETIC kin {\n", + " ~ x << (a)\n", + " ~ x -> (b)\n", + "}\n", + "\"\"\"\n", + "print(run_kinetic_visitor_and_return_derivative(ex4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 5\n", + "Reaction statements and use of `f_flux`, `b_flux` variables" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DERIVATIVE kin {\n", + " f = a*x-b*y\n", + " g = c*z\n", + " h = 0\n", + " x' = (-1*(a*x-b*y))\n", + " y' = (1*(a*x-b*y))\n", + " z' = (-1*(c*z))\n", + "}\n" + ] + } + ], + "source": [ + "ex5 = \"\"\"\n", + "STATE {\n", + " x y z\n", + "}\n", + "KINETIC kin {\n", + " ~ x <-> y (a,b)\n", + " f = f_flux - b_flux\n", + " ~ z -> (c)\n", + " g = f_flux\n", + " h = b_flux\n", + "}\n", + "\"\"\"\n", + "print(run_kinetic_visitor_and_return_derivative(ex5))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/notebooks/nmodl-linear-solver.ipynb b/docs/notebooks/nmodl-linear-solver.ipynb new file mode 100644 index 0000000000..fe6739b81f --- /dev/null +++ b/docs/notebooks/nmodl-linear-solver.ipynb @@ -0,0 +1,47 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### [wip] NMODL LINEAR solver\n", + "\n", + "`LINEAR` blocks contain a set of simultaneous equations.\n", + "\n", + "These are solved by `solve_lin_system` from [nmodl/ode.py](https://github.com/BlueBrain/nmodl/blob/master/nmodl/ode.py#L143).\n", + "\n", + "If the system is sufficiently small (by default $N\\leq3$), then Gaussian Elimination is used to directly construct the solution at compile time using SymPy to do the symbolic Gaussian Elimination.\n", + "\n", + "Otherwise at run time by LU factorization with partial pivoting." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/notebooks/nmodl-nonlinear-solver.ipynb b/docs/notebooks/nmodl-nonlinear-solver.ipynb new file mode 100644 index 0000000000..9c9de2099f --- /dev/null +++ b/docs/notebooks/nmodl-nonlinear-solver.ipynb @@ -0,0 +1,47 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### [wip] NMODL NONLINEAR solver\n", + "\n", + "`NONLINEAR` blocks contain a set of non-linear simultaneous equations.\n", + "\n", + "These are solved at runtime using Newton iteration to find the root of a function $F(X)$ with Jacobian $J = \\frac{\\partial F(X)_i}{\\partial X_j}$, where $F$ and $J$ are constructed by the [solve_non_lin_system](https://github.com/BlueBrain/nmodl/blob/master/nmodl/ode.py#L209) python routine which uses SymPy to analytically differentiate $F$ to find the Jacobian.\n", + "\n", + "The Newton solver is called `newton_solver` and is implemented in [src/solver/newton](https://github.com/BlueBrain/nmodl/blob/master/src/solver/newton/newton.hpp#L33) using the Eigen matrix algebra library.\n", + "\n", + "A fall-back solution if the analytic Jacobian is not available is to use the `newton_numerical_diff_solver` variant of this solver that uses a finite difference approximation to estimate the Jacobian numerically" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/notebooks/nmodl-odes-overview.ipynb b/docs/notebooks/nmodl-odes-overview.ipynb new file mode 100644 index 0000000000..9b49adbcf1 --- /dev/null +++ b/docs/notebooks/nmodl-odes-overview.ipynb @@ -0,0 +1,102 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### NMODL integration of ODEs\n", + "\n", + "This is an overview of\n", + " - the different ways a user can specify the equations that define the system they want to simulate in the MOD file\n", + " - how these equations can be related to each other\n", + " - how these equations are solved in NMODL\n", + "***\n", + "The user can specify information about the system in a variety of ways:\n", + " - A high level way to describe a system is to specify a Mass Action Kinetic scheme of reaction equations in a `KINETIC` block\n", + " - Alternatively a system of ODEs can be specified in a `DERIVATIVE` block (any kinetic scheme can also be written as a system of ODEs)\n", + " - Finally a system of linear/non-linear algebraic equations can be specified in a `LINEAR`/`NONLINEAR` block (a numerical integration scheme, such as backwards Euler, transforms a system of ODEs into a system of linear or non-linear equations)\n", + "\n", + "To reduce duplication of functionality for dealing with these related systems, we implement a hierarchy of transformations:\n", + " - `KINETIC` blocks of reaction statements are translated to `DERIVATIVE` blocks of equivalent ODE systems using the law of Mass Action\n", + " - `DERIVATIVE` blocks of ODEs are translated to `(NON)LINEAR` blocks of algebraic equations using a numerical integration scheme\n", + "\n", + "After these transformations we are left with only `LINEAR`/`NONLINEAR` blocks that are then solved numerically (by Gaussian Elimination, LU factorization or Newton iteration)\n", + " ***\n", + "`KINETIC` block\n", + " - Mass Action kinetics: a set of reaction equations with associated reaction rates\n", + " - converted to a `DERIVATIVE` blocking containing an equivalent system of ODEs using the law of Mass Action\n", + " - see the [nmodl-kinetic-schemes](nmodl-kinetic-schemes.ipynb) notebook for more details\n", + "***\n", + "`DERIVATIVE` block\n", + " - system of ODEs & associated solve method: `cnexp`, `sparse`, `derivimplicit` or `euler`\n", + " - `cnexp`\n", + " - applicable if ODEs are linear & independent\n", + " - exact analytic integration\n", + " - `sparse`\n", + " - applicable if ODEs are linear & coupled\n", + " - backwards Euler numerical integration scheme: $f(t+\\Delta t) = f(t) + dt f'(t+ \\Delta t)$\n", + " - results in a linear algebraic system to solve\n", + " - numerically stable\n", + " - $\\mathcal{O}(\\Delta t)$ integration error\n", + " - `derivimplcit`\n", + " - always applicable\n", + " - backwards Euler numerical integration scheme: $f(t+\\Delta t) = f(t) + dt f'(t+ \\Delta t)$\n", + " - results in a non-linear algebraic system to solve\n", + " - numerically stable\n", + " - $\\mathcal{O}(\\Delta t)$ integration error\n", + " - `euler`\n", + " - always applicable\n", + " - forwards Euler numerical integration scheme: $f(t+\\Delta t) = f(t) + \\Delta t f'(t)$\n", + " - numerically unstable\n", + " - $\\mathcal{O}(\\Delta t)$ integration error\n", + " - not recommended due to instability of scheme\n", + " - see the [nmodl-sympy-solver](nmodl-sympy-solver.ipynb) notebook for more details\n", + "***\n", + "`LINEAR` block\n", + " - system of linear algebraic equations\n", + " - for small systems ($N<=3$)\n", + " - solve at compile time by Gaussian elimination\n", + " - for larger systems\n", + " - solve at run-time by LU factorization with partial pivoting\n", + " - see the [nmodl-linear-solver](nmodl-linear-solver.ipynb) notebook for more details \n", + "***\n", + "`NONLINEAR` block\n", + " - system of non-linear algebraic equations\n", + " - solve by Newton iteration\n", + " - construct $F(X)$, with Jacobian $J(X)=\\frac{\\partial F_i}{\\partial X_j}$\n", + " - such that desired solution $X^*$ satisfies condition $F(X^*) = 0$\n", + " - iterative solution given by $X_{n+1} = X_n + J(X_n)^{-1} F(X_n)$\n", + " - see the [nmodl-nonlinear-solver](nmodl-nonlinear-solver.ipynb) notebook for more details \n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/notebooks/nmodl-python-tutorial.ipynb b/docs/notebooks/nmodl-python-tutorial.ipynb index 362dd2a282..53f790f783 100644 --- a/docs/notebooks/nmodl-python-tutorial.ipynb +++ b/docs/notebooks/nmodl-python-tutorial.ipynb @@ -15,8 +15,9 @@ "metadata": {}, "outputs": [], "source": [ - "from IPython.display import display, Javascript, HTML\n", - "import json" + "import json\n", + "\n", + "from IPython.display import HTML, Javascript, display" ] }, { @@ -259,7 +260,7 @@ } ], "source": [ - "Javascript(filename='tree.js')" + "Javascript(filename=\"tree.js\")" ] }, { @@ -313,7 +314,8 @@ } ], "source": [ - "HTML(\"\"\"\n", + "HTML(\n", + " \"\"\"\n", "\n", - "\"\"\")" + "\"\"\"\n", + ")" ] }, { @@ -464,7 +467,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we can parse any valid NMODL constructs using parsing interface. First, we have to create nmodl parser object using `nmodl::Driver` and then we can use `parse_string` method :" + "Now we can parse any valid NMODL constructs using parsing interface. First, we have to create nmodl parser object using `nmodl::NmodlDriver` and then we can use `parse_string` method :" ] }, { @@ -484,7 +487,7 @@ } ], "source": [ - "driver = nmodl.Driver()\n", + "driver = nmodl.NmodlDriver()\n", "driver.parse_string(channel)" ] }, @@ -525,8 +528,9 @@ } ], "source": [ - "print ('%.100s' % modast) # only first 100 characters\n", + "print(\"%.100s\" % modast) # only first 100 characters\n", "import json\n", + "\n", "json_data = json.loads(nmodl.to_json(modast, True))\n", "json_data_expand = json.loads(nmodl.to_json(modast, True, True))" ] @@ -540,7 +544,7 @@ "data": { "application/javascript": [ "(function(element){\n", - " require(['draw_tree'], function(draw) { draw(element.get(0), {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"SUFFIX\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"CaDynamics\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"Suffix\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"ca\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"ica\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"ReadIonVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"WriteIonVar\"}], \"name\": \"Useion\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"decay\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"depth\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}], \"name\": \"Range\"}], \"name\": \"StatementBlock\"}], \"name\": \"NeuronBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"mV\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"millivolt\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"mA\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"milliamp\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"FARADAY\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"faraday\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"coulombs\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"FactorDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"molar\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"1/liter\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"mM\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"millimolar\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"um\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"micron\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}], \"name\": \"UnitBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"0.05\"}], \"name\": \"Double\"}], \"name\": \"ParamAssign\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"decay\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"80\"}], \"name\": \"Integer\"}, {\"children\": [{\"children\": [{\"name\": \"ms\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"ParamAssign\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"depth\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"0.1\"}], \"name\": \"Double\"}, {\"children\": [{\"children\": [{\"name\": \"um\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"ParamAssign\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"0.0001\"}], \"name\": \"Double\"}, {\"children\": [{\"children\": [{\"name\": \"mM\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"ParamAssign\"}], \"name\": \"ParamBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"ica\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"mA/cm2\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"DependentDef\"}], \"name\": \"DependentBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"=\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"InitialBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"mM\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"DependentDef\"}], \"name\": \"StateBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"states\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"cnexp\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"SolveBlock\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"BreakpointBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"states\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}, {\"children\": [{\"name\": \"1\"}], \"name\": \"Integer\"}], \"name\": \"PrimeName\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"=\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"-\"}], \"name\": \"UnaryOperator\"}, {\"children\": [{\"children\": [{\"name\": \"10000\"}], \"name\": \"Double\"}], \"name\": \"ParenExpression\"}], \"name\": \"UnaryExpression\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"ica\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}, {\"children\": [{\"name\": \"/\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"2\"}], \"name\": \"Double\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"FARADAY\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"depth\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ParenExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ParenExpression\"}], \"name\": \"BinaryExpression\"}, {\"children\": [{\"name\": \"-\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"-\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ParenExpression\"}, {\"children\": [{\"name\": \"/\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"decay\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"DiffEqExpression\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"DerivativeBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"foo\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"temp\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"LocalVar\"}], \"name\": \"LocalListStatement\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"foo\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"=\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"name\": \"1\"}], \"name\": \"Double\"}, {\"children\": [{\"name\": \"+\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"FunctionBlock\"}], \"name\": \"Program\"}) });\n", + " require(['draw_tree'], function(draw) { draw(element.get(0), {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"SUFFIX\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"CaDynamics\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"Suffix\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"ca\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"ica\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"ReadIonVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"WriteIonVar\"}], \"name\": \"Useion\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"decay\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"depth\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"RangeVar\"}], \"name\": \"Range\"}], \"name\": \"StatementBlock\"}], \"name\": \"NeuronBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"mV\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"millivolt\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"mA\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"milliamp\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"FARADAY\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"faraday\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"coulombs\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"FactorDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"molar\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"1/liter\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"mM\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"millimolar\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"um\"}], \"name\": \"String\"}], \"name\": \"Unit\"}, {\"children\": [{\"children\": [{\"name\": \"micron\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"UnitDef\"}], \"name\": \"UnitBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"0.05\"}], \"name\": \"Double\"}], \"name\": \"ParamAssign\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"decay\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"80\"}], \"name\": \"Integer\"}, {\"children\": [{\"children\": [{\"name\": \"ms\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"ParamAssign\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"depth\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"0.1\"}], \"name\": \"Double\"}, {\"children\": [{\"children\": [{\"name\": \"um\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"ParamAssign\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"name\": \"0.0001\"}], \"name\": \"Double\"}, {\"children\": [{\"children\": [{\"name\": \"mM\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"ParamAssign\"}], \"name\": \"ParamBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"ica\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"mA/cm2\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"DependentDef\"}], \"name\": \"DependentBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"=\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"InitialBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"mM\"}], \"name\": \"String\"}], \"name\": \"Unit\"}], \"name\": \"DependentDef\"}], \"name\": \"StateBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"states\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"name\": \"cnexp\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"SolveBlock\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"BreakpointBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"states\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}, {\"children\": [{\"name\": \"1\"}], \"name\": \"Integer\"}], \"name\": \"PrimeName\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"=\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"-\"}], \"name\": \"UnaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"10000\"}], \"name\": \"Double\"}], \"name\": \"ParenExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"UnaryExpression\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"ica\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}, {\"children\": [{\"name\": \"/\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"2\"}], \"name\": \"Double\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"FARADAY\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}, {\"children\": [{\"name\": \"*\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"depth\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"ParenExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"ParenExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}, {\"children\": [{\"name\": \"-\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"cai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"-\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"minCai\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"ParenExpression\"}], \"name\": \"WrappedExpression\"}, {\"children\": [{\"name\": \"/\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"decay\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"DiffEqExpression\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"DerivativeBlock\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"foo\"}], \"name\": \"String\"}], \"name\": \"Name\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"temp\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"LocalVar\"}], \"name\": \"LocalListStatement\"}, {\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"children\": [{\"name\": \"foo\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}, {\"children\": [{\"name\": \"=\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"1\"}], \"name\": \"Double\"}, {\"children\": [{\"name\": \"+\"}], \"name\": \"BinaryOperator\"}, {\"children\": [{\"children\": [{\"children\": [{\"name\": \"gamma\"}], \"name\": \"String\"}], \"name\": \"Name\"}], \"name\": \"VarName\"}], \"name\": \"BinaryExpression\"}], \"name\": \"WrappedExpression\"}], \"name\": \"BinaryExpression\"}], \"name\": \"ExpressionStatement\"}], \"name\": \"StatementBlock\"}], \"name\": \"FunctionBlock\"}], \"name\": \"Program\"}) });\n", " })(element);" ], "text/plain": [ @@ -553,9 +557,12 @@ } ], "source": [ - "Javascript(\"\"\"(function(element){\n", + "Javascript(\n", + " \"\"\"(function(element){\n", " require(['draw_tree'], function(draw) { draw(element.get(0), %s) });\n", - " })(element);\"\"\" % json.dumps(json_data_expand) )" + " })(element);\"\"\"\n", + " % json.dumps(json_data_expand)\n", + ")" ] }, { @@ -582,8 +589,8 @@ "metadata": {}, "outputs": [], "source": [ - "from nmodl.dsl import visitor\n", - "from nmodl.dsl import ast\n", + "from nmodl.dsl import ast, visitor\n", + "\n", "lookup_visitor = visitor.AstLookupVisitor()" ] }, @@ -612,7 +619,7 @@ "source": [ "states = lookup_visitor.lookup(modast, ast.AstNodeType.STATE_BLOCK)\n", "for state in states:\n", - " print (nmodl.to_nmodl(state))" + " print(nmodl.to_nmodl(state))" ] }, { @@ -652,30 +659,30 @@ "if odes:\n", " print(\"%d differential equation(s) exist : \" % len(odes))\n", " for ode in odes:\n", - " print (\"\\t %s \" % nmodl.to_nmodl(ode))\n", - " \n", + " print(\"\\t %s \" % nmodl.to_nmodl(ode))\n", + "\n", "if primes:\n", - " print(\"%d prime variables exist :\" % len(primes), end='')\n", + " print(\"%d prime variables exist :\" % len(primes), end=\"\")\n", " for prime in primes:\n", - " print (\" %s\" % nmodl.to_nmodl(prime), end='')\n", + " print(\" %s\" % nmodl.to_nmodl(prime), end=\"\")\n", " print()\n", "\n", "if range_vars:\n", - " print(\"%d range variables exist :\" % len(range_vars), end='')\n", + " print(\"%d range variables exist :\" % len(range_vars), end=\"\")\n", " for range_var in range_vars:\n", - " print (\" %s\" % nmodl.to_nmodl(range_var), end='')\n", + " print(\" %s\" % nmodl.to_nmodl(range_var), end=\"\")\n", " print()\n", "\n", "if parameters:\n", - " print(\"%d parameter variables exist :\" % len(parameters), end='')\n", + " print(\"%d parameter variables exist :\" % len(parameters), end=\"\")\n", " for range_var in range_vars:\n", - " print (\" %s\" % nmodl.to_nmodl(range_var), end='')\n", + " print(\" %s\" % nmodl.to_nmodl(range_var), end=\"\")\n", " print()\n", - " \n", + "\n", "if units:\n", - " print(\"%d units uses :\" % len(units), end='')\n", + " print(\"%d units uses :\" % len(units), end=\"\")\n", " for unit in units:\n", - " print (\" %s\" % nmodl.to_nmodl(unit), end='')" + " print(\" %s\" % nmodl.to_nmodl(unit), end=\"\")" ] }, { @@ -700,7 +707,7 @@ ], "source": [ "functions = lookup_visitor.lookup(modast, ast.AstNodeType.FUNCTION_BLOCK)\n", - "function = functions[0] # first function\n", + "function = functions[0] # first function\n", "\n", "# expression statements include assignments\n", "new_lookup_visitor = visitor.AstLookupVisitor(ast.AstNodeType.EXPRESSION_STATEMENT)\n", @@ -710,7 +717,7 @@ "statements = new_lookup_visitor.get_nodes()\n", "\n", "for statement in statements:\n", - " print (nmodl.to_nmodl(statement))" + " print(nmodl.to_nmodl(statement))" ] }, { @@ -795,7 +802,7 @@ "table = modast.get_symbol_table()\n", "table_s = str(table)\n", "\n", - "print (table_s)" + "print(table_s)" ] }, { @@ -819,8 +826,8 @@ } ], "source": [ - "cai = table.lookup('cai')\n", - "print (cai)" + "cai = table.lookup(\"cai\")\n", + "print(cai)" ] }, { @@ -855,7 +862,7 @@ "source": [ "range_vars = table.get_variables_with_properties(symtab.NmodlType.range_var)\n", "for var in range_vars:\n", - " print (var)" + " print(var)" ] }, { @@ -880,9 +887,11 @@ } ], "source": [ - "ions_var = table.get_variables_with_properties(symtab.NmodlType.read_ion_var | symtab.NmodlType.write_ion_var, False)\n", + "ions_var = table.get_variables_with_properties(\n", + " symtab.NmodlType.read_ion_var | symtab.NmodlType.write_ion_var, False\n", + ")\n", "for var in ions_var:\n", - " print (var)" + " print(var)" ] }, { @@ -915,11 +924,12 @@ "source": [ "from nmodl.dsl import ast, visitor\n", "\n", - "class DoubleVisitor(visitor.AstVisitor):\n", "\n", + "class DoubleVisitor(visitor.AstVisitor):\n", " def visit_double(self, node):\n", - " print (node.eval()) # or, can use nmodl.to_nmodl(node) \n", - " \n", + " print(node.eval()) # or, can use nmodl.to_nmodl(node)\n", + "\n", + "\n", "d_visitor = DoubleVisitor()\n", "modast.accept(d_visitor)" ] @@ -955,11 +965,10 @@ ], "source": [ "class ParameterVisitor(visitor.AstVisitor):\n", - " \n", " def __init__(self):\n", " visitor.AstVisitor.__init__(self)\n", " self.in_parameter = False\n", - " \n", + "\n", " def visit_param_block(self, node):\n", " self.in_parameter = True\n", " node.visit_children(self)\n", @@ -967,15 +976,16 @@ "\n", " def visit_name(self, node):\n", " if self.in_parameter:\n", - " print (nmodl.to_nmodl(node))\n", - " \n", + " print(nmodl.to_nmodl(node))\n", + "\n", " def visit_double(self, node):\n", " if self.in_parameter:\n", - " print (node.eval())\n", + " print(node.eval())\n", "\n", " def visit_integer(self, node):\n", " if self.in_parameter:\n", - " print (node.eval())\n", + " print(node.eval())\n", + "\n", "\n", "param_visitor = ParameterVisitor()\n", "modast.accept(param_visitor)" @@ -998,7 +1008,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.7.2" } }, "nbformat": 4, diff --git a/docs/notebooks/nmodl-sympy-conductance.ipynb b/docs/notebooks/nmodl-sympy-conductance.ipynb new file mode 100644 index 0000000000..f205a3584d --- /dev/null +++ b/docs/notebooks/nmodl-sympy-conductance.ipynb @@ -0,0 +1,455 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### NMODL CONDUCTANCE\n", + "\n", + "This notebook described the `CONDUCTANCE` keyword in NEURON, how it is implemented in NMODL, and shows some examples of the output generated by NMODL in different situations.\n", + "\n", + "For a more general tutorial on using the NMODL python interface, please see the [tutorial notebook](nmodl-python-tutorial.ipynb)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Introduction\n", + "\n", + "Motivation:\n", + " - during a NEURON simulation, a number of currents $I$ may be generated\n", + " - at each time step, for each current $I$ the corresponding conductance $dI/dV$ needs to be calculated\n", + " - by default in NEURON this is approximated by a forwards difference: $dI/dV \\simeq (I(v+\\Delta v)-I(v))/\\Delta v$, with $\\Delta v = 0.001$\n", + " - this introduces an $\\mathcal{O}(\\Delta v)$ numerical error\n", + " - it also requires two current calculations, which may be computationally inefficient\n", + "\n", + "Solution:\n", + " - the `CONDUCTANCE` keyword was added to the NMODL language\n", + " - this allows the user to manually specify the analytic expression for the conductance in the MOD file\n", + " - during the simulation, instead of the numerical differentiation, the user supplied expression is used\n", + " - this solves the problem, but requires additional effort from the user\n", + " - it also opens up room for user error: an incorrect expression will still run but the results will not be correct\n", + "\n", + "SymPy improvement:\n", + " - the currents in the mod file are differentiated analytically using SymPy\n", + " - the corresponding `CONDUCTANCE` statements are generated automatically\n", + " - no additional input required from the user\n", + " - avoids the possibility of user error" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Implementation\n", + "\n", + "The `SympyConductanceVisitor` is defined in [src/visitors/sympy_conductance_visitor.hpp](https://github.com/BlueBrain/nmodl/blob/master/src/visitors/sympy_conductance_visitor.hpp), and it makes use of the python function [differentiate2c](https://github.com/BlueBrain/nmodl/blob/master/nmodl/ode.py#L377) to perform the analytic differentation using the [SymPy](https://www.sympy.org/en/index.html) symbolic math Python library.\n", + "\n", + "* For each ion write statement $i = \\dots$ in the BREAKPOINT block\n", + " * Differentiate to find the conductance $g_i=di/dv$\n", + " * If this $g_i$ coincides with an existing variable, e.g. $g$, add to BREAKPOINT the statement:\n", + " * CONDUCTANCE g USEION ion_name\n", + " * If not, also need to declare and asign a variable for the calculated conductance:\n", + " * LOCAL g_i_0\n", + " * CONDUCTANCE g_i_0 USEION ion_name\n", + " * g_i_0 = ...\n", + " * But if there is an existing CONDUCTANCE statement, then do not modify it\n", + "\n", + "\n", + "* It may be the case that a variable in the write statement $i = \\dots$ itself depends on $v$, so to take this into account:\n", + " * an inlining visitor is first ran, after which all variable assignments occur within the BREAKPOINT block\n", + " * each preceeding expression is analysed in reverse order for $v$ dependence \n", + " * if it depends on $v$, the rhs of the expression is substituted for the lhs in all following statements\n", + " * the end result is a (complicated) expression $i = ...$ where all v dependence is explicit\n", + " * this is then differentiated w.r.t $v$ to give the conductance\n", + " * it then checks if this expression is equivalent to an existing variable\n", + " * for this step it is necessary to also substitute all non-$v$-dependent expressions on both sides & simplify\n", + "\n", + "#### Implementation Tests\n", + "\n", + " - The unit tests may be helpful to understand what these functions are doing\n", + " - `SympyConductanceVisitor` tests are located in [test/visitor/visitor.cpp](https://github.com/BlueBrain/nmodl/blob/master/test/visitor/visitor.cpp#L3261)\n", + " - `differentiate2c` tests are located in [test/ode/test_ode.py](https://github.com/BlueBrain/nmodl/blob/master/test/ode/test_ode.py#L56)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Examples\n", + "Here are some examples of generated CONDUCTANCE statements for a variety of sample mod files." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import nmodl.dsl as nmodl\n", + "\n", + "\n", + "def run_conductance_visitor_and_return_breakpoint(mod_string):\n", + " # parse NMDOL file (supplied as a string) into AST\n", + " driver = nmodl.NmodlDriver()\n", + " driver.parse_string(mod_string)\n", + " AST = driver.ast()\n", + " # run SymtabVisitor to generate Symbol Table\n", + " nmodl.symtab.SymtabVisitor().visit_program(AST)\n", + " # constant folding, inlining & local variable renaming passes\n", + " nmodl.visitor.ConstantFolderVisitor().visit_program(AST)\n", + " nmodl.visitor.InlineVisitor().visit_program(AST)\n", + " nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)\n", + "\n", + " # run CONDUCTANCE visitor\n", + " nmodl.visitor.SympyConductanceVisitor().visit_program(AST)\n", + " # return new BREAKPOINT block\n", + " return nmodl.to_nmodl(\n", + " nmodl.visitor.AstLookupVisitor().lookup(\n", + " AST, nmodl.ast.AstNodeType.BREAKPOINT_BLOCK\n", + " )[0]\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 1\n", + " - simple USEION statement, conductance equal to existing variable\n", + " - add CONDUCTANCE statement using existing variable" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " CONDUCTANCE gna USEION na\n", + " ina = gna*(v-ena)\n", + "}\n" + ] + } + ], + "source": [ + "ex1 = \"\"\"\n", + "NEURON {\n", + " USEION na READ ena WRITE ina\n", + " RANGE gna\n", + "}\n", + "BREAKPOINT {\n", + " ina = gna*(v - ena)\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 2\n", + " - simple USEION statement, conductance not equal to existing variable\n", + " - declare new local variable\n", + " - assign conductance to it\n", + " - add CONDUCTANCE statement" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " LOCAL g_na_0\n", + " CONDUCTANCE g_na_0 USEION na\n", + " g_na_0 = 0.1*gna\n", + " ina = 0.1*gna*(v-ena)\n", + "}\n" + ] + } + ], + "source": [ + "ex2 = \"\"\"\n", + "NEURON {\n", + " USEION na READ ena WRITE ina\n", + " RANGE gna\n", + "}\n", + "BREAKPOINT {\n", + " ina = 0.1*gna*(v - ena)\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 3\n", + " - simple NONSPECIFIC_CURRENT statement, conductance equal to existing variable\n", + " - add CONDUCTANCE statement using existing variable" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " CONDUCTANCE g\n", + " i = g*v\n", + "}\n" + ] + } + ], + "source": [ + "ex3 = \"\"\"\n", + "NEURON {\n", + " NONSPECIFIC_CURRENT i\n", + " RANGE g\n", + "}\n", + "BREAKPOINT {\n", + " i = g*v\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 4\n", + " - non-linear NONSPECIFIC_CURRENT statement, conductance not equal to existing variable\n", + " - declare new local variable\n", + " - assign conductance to it\n", + " - add CONDUCTANCE statement" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " LOCAL g__0\n", + " CONDUCTANCE g__0\n", + " g__0 = g+2*v\n", + " i = g*v+v*v\n", + "}\n" + ] + } + ], + "source": [ + "ex4 = \"\"\"\n", + "NEURON {\n", + " NONSPECIFIC_CURRENT i\n", + " RANGE g\n", + "}\n", + "BREAKPOINT {\n", + " i = g*v + v*v\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 5\n", + " - several current statements, conductance equal to existing variables\n", + " - add CONDUCTANCE statements" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " CONDUCTANCE gl\n", + " CONDUCTANCE gk USEION k\n", + " CONDUCTANCE gna USEION na\n", + " gna = gnabar*m*m*m*h\n", + " ina = gna*(v-ena)\n", + " gk = gkbar*n*n*n*n\n", + " ik = gk*(v-ek)\n", + " il = gl*(v-el)\n", + "}\n" + ] + } + ], + "source": [ + "ex5 = \"\"\"\n", + "NEURON {\n", + " USEION na READ ena WRITE ina\n", + " USEION k READ ek WRITE ik\n", + " NONSPECIFIC_CURRENT il\n", + " RANGE gnabar, gkbar, gl, el, gna, gk\n", + "}\n", + "STATE {\n", + " m n h\n", + "}\n", + "BREAKPOINT {\n", + " gna = gnabar*m*m*m*h\n", + " ina = gna*(v - ena)\n", + " gk = gkbar*n*n*n*n\n", + " ik = gk*(v - ek)\n", + " il = gl*(v - el)\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex5))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 6\n", + " - current contains variables that depend on $v$, conductance equal to existing variable `x3`\n", + " - substitute all variables with $v$-dependence, differentiate to find conductance\n", + " - compare result to each existing variable (after substituting all preceeding declarations on both sides)\n", + " - identify that expression for conductance is equivalent to `x3`\n", + " - add CONDUCTANCE statement using this existing variable" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " CONDUCTANCE x3 USEION na\n", + " x1 = 0.2+3*v\n", + " x2 = v*v\n", + " x3 = 3*v*v+5*v-1.3\n", + " gna = x1+x2\n", + " ina = gna*(v-0.5)\n", + "}\n" + ] + } + ], + "source": [ + "ex6 = \"\"\"\n", + "NEURON {\n", + " USEION na READ ena WRITE ina\n", + " RANGE gna, x1, x2, x3\n", + "}\n", + "BREAKPOINT {\n", + " x1 = 0.2+3*v\n", + " x2 = v*v\n", + " x3 = 3*v*v+5*v-1.3\n", + " gna = x1 + x2\n", + " ina = gna*(v-0.5)\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex6))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Ex. 7\n", + " - current contains variables that depend on $v$, conductance not equal to existing variable\n", + " - substitute all variables with $v$-dependence, differentiate to find conductance\n", + " - compare result to each existing variable, no equivalent expression found\n", + " - declare new local variable, assign conductance to it\n", + " - add CONDUCTANCE statement using this new variable" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BREAKPOINT {\n", + " LOCAL g_na_0\n", + " CONDUCTANCE g_na_0 USEION na\n", + " g_na_0 = 3*pow(v, 2)+5*v-1.3\n", + " x1 = 0.2+3*v\n", + " x2 = v*v\n", + " gna = x1+x2\n", + " ina = gna*(v-0.5)\n", + "}\n" + ] + } + ], + "source": [ + "ex7 = \"\"\"\n", + "NEURON {\n", + " USEION na READ ena WRITE ina\n", + " RANGE gna, x1, x2\n", + "}\n", + "BREAKPOINT {\n", + " x1 = 0.2+3*v\n", + " x2 = v*v\n", + " gna = x1 + x2\n", + " ina = gna*(v-0.5)\n", + "}\n", + "\"\"\"\n", + "print(run_conductance_visitor_and_return_breakpoint(ex7))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/notebooks/nmodl-python-sympy-examples.ipynb b/docs/notebooks/nmodl-sympy-solver.ipynb similarity index 66% rename from docs/notebooks/nmodl-python-sympy-examples.ipynb rename to docs/notebooks/nmodl-sympy-solver.ipynb index 14980bbbbc..ac4ebfcccf 100644 --- a/docs/notebooks/nmodl-python-sympy-examples.ipynb +++ b/docs/notebooks/nmodl-sympy-solver.ipynb @@ -4,13 +4,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# NMODL SymPy Visitor Examples\n", + "### [wip] NMODL SympySolver\n", "\n", - "Some examples of the use of SymPy within NMODL to\n", - "- solve differential equations to generate solutions for the CNEXP solver (`SympySolverVisitor`)\n", - "- differentiate current expressions to generate CONDUCTANCE statements (`SympyConductanceVisitor`)\n", + "This notebook describes the implementation of the `SympySolverVisitor`, which solves the systems of ODEs defined in `DERIVATIVE` blocks.\n", "\n", - "Please see the [tutorial notebook](nmodl-python-tutorial.ipynb) for a more general tutorial on using the NMODL python interface, including installation instructions." + "For a higher level overview of the approach to solving ODEs in NMODL, please see the [nmodl-odes-overview](nmodl-odes-overview.ipynb) notebook. \n", + "\n", + "For a more general tutorial on using the NMODL python interface, please see the [tutorial notebook](nmodl-python-tutorial.ipynb)." ] }, { @@ -138,7 +138,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## ODE Solver example" + "### ODE Solver example" ] }, { @@ -149,7 +149,7 @@ "source": [ "def parse_mod_to_ast(mod_string):\n", " # parse NMDOL file (supplied as a string)\n", - " driver = nmodl.Driver()\n", + " driver = nmodl.NmodlDriver()\n", " driver.parse_string(mod_string)\n", " modast = driver.ast()\n", " # run SymtabVisitor to generate Symbol Table\n", @@ -158,6 +158,7 @@ " # return AST\n", " return modast\n", "\n", + "\n", "lookup_visitor = visitor.AstLookupVisitor()" ] }, @@ -253,9 +254,9 @@ "text": [ "DERIVATIVE states {\n", " rates(v, celsius)\n", - " m = minf+(m-minf)*exp(-dt/mtau)\n", - " h = hinf+(h-hinf)*exp(-dt/htau)\n", - " n = ninf+(n-ninf)*exp(-dt/ntau)\n", + " m = minf-(-m+minf)*exp(-dt/mtau)\n", + " h = hinf-(-h+hinf)*exp(-dt/htau)\n", + " n = ninf-(-n+ninf)*exp(-dt/ntau)\n", "}\n" ] } @@ -311,131 +312,6 @@ " nmodl.to_nmodl(lookup_visitor.lookup(modast, ast.AstNodeType.DERIVATIVE_BLOCK)[0])\n", ")" ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## CONDUCTANCE example" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The CONDUCTANCE keyword has been introduced to NEURON as well as CoreNEURON. If the i/v relation is ohmic in BREAKPOINT block then one can use CONDUCTANCE keyword for efficiency." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "USEION na READ ena WRITE ina\n", - "USEION k READ ek WRITE ik\n", - "NONSPECIFIC_CURRENT il\n" - ] - } - ], - "source": [ - "modast = parse_mod_to_ast(channel)\n", - "# print USEION and NONSPECIFIC current statements\n", - "for node in lookup_visitor.lookup(modast, ast.AstNodeType.USEION):\n", - " print(nmodl.to_nmodl(node))\n", - "for node in lookup_visitor.lookup(modast, ast.AstNodeType.NONSPECIFIC):\n", - " print(nmodl.to_nmodl(node))" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "BREAKPOINT {\n", - " SOLVE states METHOD cnexp\n", - " gna = gnabar*m*m*m*h\n", - " ina = gna*(v-ena)\n", - " gk = gkbar*n*n*n*n\n", - " ik = gk*(v-ek)\n", - " il = gl*(v-el)\n", - "}\n" - ] - } - ], - "source": [ - "# print BREAKPOINT\n", - "print(\n", - " nmodl.to_nmodl(lookup_visitor.lookup(modast, ast.AstNodeType.BREAKPOINT_BLOCK)[0])\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If we run `SympyConductanceVisitor`, it does the following:\n", - "\n", - "* For each ion write statement $i = \\dots$ in the BREAKPOINT block\n", - " * Differentiate to find the conductance $g_i=di/dv$\n", - " * If this $g_i$ coincides with an existing variable, e.g. $g$, add to BREAKPOINT the statement:\n", - " * CONDUCTANCE g USEION ion_name\n", - " * If not, also need to declare and asign a variable for the calculated conductance:\n", - " * LOCAL g_i_0\n", - " * CONDUCTANCE g_i_0 USEION ion_name\n", - " * g_i_0 = ...\n", - " * But if there is an existing CONDUCTANCE statement, then do not modify it\n", - "\n", - "\n", - "* NOTE: currently we just differentiate the equation, assuming that the variables do not depend on v\n", - " * Ideally we should do something like:\n", - " * check after inlining which variables are written to in BREAKPOINT block\n", - " * get rhs of each such expression, and substitute it (recursively) for the lhs in SymPy\n", - " * this should give a (complicated) expression $i = ...$ where all v dependence is explicit\n", - " * then differentiate this w.r.t v\n", - " * then simplify and try to write the result in terms of an existing variable\n", - "\n", - "Here is the BREAKPOINT block after running `SympyConductanceVisitor`:" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "BREAKPOINT {\n", - " CONDUCTANCE gna USEION na\n", - " CONDUCTANCE gl\n", - " CONDUCTANCE gk USEION k\n", - " SOLVE states METHOD cnexp\n", - " gna = gnabar*m*m*m*h\n", - " ina = gna*(v-ena)\n", - " gk = gkbar*n*n*n*n\n", - " ik = gk*(v-ek)\n", - " il = gl*(v-el)\n", - "}\n" - ] - } - ], - "source": [ - "sympy_conductance_visitor = visitor.SympyConductanceVisitor()\n", - "sympy_conductance_visitor.visit_program(modast)\n", - "# print BREAKPOINT block\n", - "print(\n", - " nmodl.to_nmodl(lookup_visitor.lookup(modast, ast.AstNodeType.BREAKPOINT_BLOCK)[0])\n", - ")" - ] } ], "metadata": { @@ -454,7 +330,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.7.2" } }, "nbformat": 4, diff --git a/src/language/parser.py b/src/language/parser.py index 5120cb6177..249e6fcc7b 100644 --- a/src/language/parser.py +++ b/src/language/parser.py @@ -191,7 +191,7 @@ def parse_file(self): with open(self.filename, 'r') as stream: try: - rules = yaml.load(stream) + rules = yaml.safe_load(stream) _, nodes = self.parse_yaml_rules(rules) except yaml.YAMLError as e: print("Error while parsing YAML definition file {0} : {1}".format( diff --git a/src/language/templates/pyvisitor.cpp b/src/language/templates/pyvisitor.cpp index 993f7f4958..8126d53083 100644 --- a/src/language/templates/pyvisitor.cpp +++ b/src/language/templates/pyvisitor.cpp @@ -12,6 +12,10 @@ #include "pybind/pybind_utils.hpp" #include "pybind/pyvisitor.hpp" +#include "visitors/constant_folder_visitor.hpp" +#include "visitors/inline_visitor.hpp" +#include "visitors/kinetic_block_visitor.hpp" +#include "visitors/local_var_rename_visitor.hpp" #include "visitors/lookup_visitor.hpp" #include "visitors/nmodl_visitor.hpp" #include "visitors/sympy_conductance_visitor.hpp" @@ -43,15 +47,36 @@ static const char *nmodlprintvisitor_class = R"( NmodlPrintVisitor class )"; -static const char *sympysolvervisitor_class = R"( +static const char *constantfoldervisitor_class = R"( - SympySolverVisitor class + ConstantFolderVisitor class +)"; + +static const char *inlinevisitor_class = R"( + + InlineVisitor class +)"; + +static const char *kineticblockvisitor_class = R"( + + KineticBlockVisitor class +)"; + +static const char *localvarrenamevisitor_class = R"( + + LocalVarRenameVisitor class )"; static const char *sympyconductancevisitor_class = R"( SympyConductanceVisitor class )"; + +static const char *sympysolvervisitor_class = R"( + + SympySolverVisitor class +)"; + } #pragma clang diagnostic push @@ -137,13 +162,29 @@ void init_visitor_module(py::module& m) { {% if loop.last -%};{% endif %} {% endfor %} - py::class_ sympy_solver_visitor(m_visitor, "SympySolverVisitor", docstring::sympysolvervisitor_class); - sympy_solver_visitor.def(py::init(), py::arg("use_pade_approx")=false) - .def("visit_program", &SympySolverVisitor::visit_program); + py::class_ constant_folder_visitor(m_visitor, "ConstantFolderVisitor", docstring::constantfoldervisitor_class); + constant_folder_visitor.def(py::init<>()) + .def("visit_program", &ConstantFolderVisitor::visit_program); + + py::class_ inline_visitor(m_visitor, "InlineVisitor", docstring::inlinevisitor_class); + inline_visitor.def(py::init<>()) + .def("visit_program", &InlineVisitor::visit_program); + + py::class_ kinetic_block_visitor(m_visitor, "KineticBlockVisitor", docstring::kineticblockvisitor_class); + kinetic_block_visitor.def(py::init<>()) + .def("visit_program", &KineticBlockVisitor::visit_program); + + py::class_ local_var_rename_visitor(m_visitor, "LocalVarRenameVisitor", docstring::localvarrenamevisitor_class); + local_var_rename_visitor.def(py::init<>()) + .def("visit_program", &LocalVarRenameVisitor::visit_program); py::class_ sympy_conductance_visitor(m_visitor, "SympyConductanceVisitor", docstring::sympyconductancevisitor_class); sympy_conductance_visitor.def(py::init<>()) .def("visit_program", &SympyConductanceVisitor::visit_program); + + py::class_ sympy_solver_visitor(m_visitor, "SympySolverVisitor", docstring::sympysolvervisitor_class); + sympy_solver_visitor.def(py::init(), py::arg("use_pade_approx")=false) + .def("visit_program", &SympySolverVisitor::visit_program); } #pragma clang diagnostic pop