Merge pull request #247 from DimitriPlotnikov/master

Enable more usefull appliaction of the ODE aliases

models/iaf_psc_exp.nestml

@@ -1,3 +1,79 @@
+/*
+  Name: iaf_psc_exp_nestml - Leaky integrate-and-fire neuron model with exponential
+                     PSCs.
+
+  Description:
+  iaf_psc_exp_nestml is an implementation of a leaky integrate-and-fire model
+  with exponential shaped postsynaptic currents (PSCs) according to [1].
+  Thus, postsynaptic currents have an infinitely short rise time.
+
+  The threshold crossing is followed by an absolute refractory period (t_ref)
+  during which the membrane potential is clamped to the resting potential
+  and spiking is prohibited.
+
+  The linear subthresold dynamics is integrated by the Exact
+  Integration scheme [2]. The neuron dynamics is solved on the time
+  grid given by the computation step size. Incoming as well as emitted
+  spikes are forced to that grid.
+
+  An additional state variable and the corresponding differential
+  equation represents a piecewise constant external current.
+
+  The general framework for the consistent formulation of systems with
+  neuron like dynamics interacting by point events is described in
+  [2]. A flow chart can be found in [3].
+
+  Remarks:
+  The present implementation uses individual variables for the
+  components of the state vector and the non-zero matrix elements of
+  the propagator.  Because the propagator is a lower triangular matrix
+  no full matrix multiplication needs to be carried out and the
+  computation can be done "in place" i.e. no temporary state vector
+  object is required.
+
+  The template support of recent C++ compilers enables a more succinct
+  formulation without loss of runtime performance already at minimal
+  optimization levels. A future version of iaf_psc_exp_nestml will probably
+  address the problem of efficient usage of appropriate vector and
+  matrix objects.
+
+  Remarks: If tau_m is very close to tau_syn_ex or tau_syn_in, the model
+  will numerically behave as if tau_m is equal to tau_syn_ex or
+  tau_syn_in, respectively, to avoid numerical instabilities.
+  For details, please see IAF_Neruons_Singularity.ipynb in the
+  NEST source code (docs/model_details).
+
+  iaf_psc_exp_nestml can handle current input in two ways: Current input
+  through receptor_type 0 are handled as stepwise constant current
+  input as in other iaf models, i.e., this current directly enters
+  the membrane potential equation. Current input through
+  receptor_type 1, in contrast, is filtered through an exponential
+  kernel with the time constant of the excitatory synapse,
+  tau_syn_ex. For an example application, see [4].
+
+  References:
+  [1] Misha Tsodyks, Asher Uziel, and Henry Markram (2000) Synchrony Generation
+  in Recurrent Networks with Frequency-Dependent Synapses, The Journal of
+  Neuroscience, 2000, Vol. 20 RC50 p. 1-5
+  [2] Rotter S & Diesmann M (1999) Exact simulation of time-invariant linear
+  systems with applications to neuronal modeling. Biologial Cybernetics
+  81:381-402.
+  [3] Diesmann M, Gewaltig M-O, Rotter S, & Aertsen A (2001) State space
+  analysis of synchronous spiking in cortical neural networks.
+  Neurocomputing 38-40:565-571.
+  [4] Schuecker J, Diesmann M, Helias M (2015) Modulated escape from a
+  metastable state driven by colored noise.
+  Physical Review E 92:052119
+
+  Sends: SpikeEvent
+
+  Receives: SpikeEvent, CurrentEvent, DataLoggingRequest
+
+  SeeAlso: iaf_psc_exp_ps
+
+  FirstVersion: March 2006
+  Author: Moritz Helias
+*/
 neuron iaf_psc_exp_nestml:
 
   state:
@@ -8,19 +84,21 @@ neuron iaf_psc_exp_nestml:
   equations:
     shape I_shape_in = exp(-1/tau_syn_in*t)
     shape I_shape_ex = exp(-1/tau_syn_ex*t)
-    V_abs' = -1/Tau * V_abs + 1/C_m * (I_sum(I_shape_in, in_spikes) + I_sum(I_shape_ex, ex_spikes) + I_e + currents)
+    I_syn mV = (I_sum(I_shape_in, in_spikes) + I_sum(I_shape_ex, ex_spikes) + I_e + currents)
+    V_abs' = -1/tau_m * V_abs + 1/C_m * I_syn
   end
 
   parameter:
-    C_m     pF = 250   # Capacity of the membrane
-    Tau     ms = 10    # Membrane time constant.
-    tau_syn_in ms = 2     # Time constant of synaptic current.
-    tau_syn_ex ms = 2     # Time constant of synaptic current.
-    t_ref   ms = 2     # Refractory period.
-    E_L     mV = -70  # Resting potential.
-    alias V_reset mV = -70 - E_L
-    alias Theta   mV = -55 - E_L
-    I_e     pA = 0     # External current.
+    C_m pF = 250pF       # Capacity of the membrane
+    tau_m ms = 10ms      # Membrane time constant.
+    tau_syn_in ms = 2ms  # Time constant of synaptic current.
+    tau_syn_ex ms = 2ms  # Time constant of synaptic current.
+    t_ref ms = 2ms       # Duration of refractory period
+    E_L  mV = -70mV      # Resting potential.
+    alias V_reset mV = -70mV - E_L # reset value of the membrane potential
+    alias Theta   mV = -55mV - E_L # Threshold, RELATIVE TO RESTING POTENTAIL(!).
+                                   # I.e. the real threshold is (E_L_+V_th_)
+    I_e     pA = 0pA     # External current.
   end
 
   internal:
@@ -37,10 +115,10 @@ neuron iaf_psc_exp_nestml:
   output: spike
 
   update:
-    if r == 0: # not refractory
+    if r == 0: # neuron not refractory, so evolve V
         integrate(V_abs)
     else:
-        r = r - 1
+        r = r - 1 # neuron is absolute refractory
     end
 
     if V_abs >= Theta: # threshold crossing

src/main/java/org/nest/codegeneration/sympy/ODETransformer.java

@@ -21,17 +21,17 @@
  */
 public class ODETransformer {
   // this function is used in freemarker templates und must be public
-  public static ASTEquation replaceSumCalls(final ASTEquation astOde) {
+  public static <T extends ASTNode> T replaceSumCalls(final T astOde) {
     // since the transformation replaces the call inplace, make a copy to preserve the information for further steps
-    final ASTEquation workingCopy = astOde.deepClone();
+    final T workingCopy = (T) astOde.deepClone(); // IT is OK, since the deepClone returns T
     final List<ASTFunctionCall> functions = get_sumFunctionCalls(workingCopy);
 
     functions.forEach(node -> replaceFunctionCallThroughFirstArgument(workingCopy, node));
     return workingCopy;
   }
 
 
-  public static List<ASTFunctionCall> get_sumFunctionCalls(final ASTNode workingCopy) {
+  private static List<ASTFunctionCall> get_sumFunctionCalls(final ASTNode workingCopy) {
     return ASTUtils.getAll(workingCopy, ASTFunctionCall.class)
           .stream()
           .filter(astFunctionCall ->
@@ -40,14 +40,14 @@ public static List<ASTFunctionCall> get_sumFunctionCalls(final ASTNode workingCo
           .collect(Collectors.toList());
   }
 
-  public static ASTExpr replaceSumCalls(final ASTExpr astExpr) {
+  /*public static ASTExpr replaceSumCalls(final ASTExpr astExpr) {
     // since the transformation replaces the call inplace, make a copy to preserve the information for further steps
     final ASTExpr workingCopy = astExpr.deepClone();
     final List<ASTFunctionCall> functions = get_sumFunctionCalls(workingCopy);
 
     functions.forEach(node -> replaceFunctionCallThroughFirstArgument(workingCopy, node));
     return workingCopy;
-  }
+  }*/
 
   public static List<ASTFunctionCall> getCondSumFunctionCall(final ASTNode workingCopy) {
     return ASTUtils.getAll(workingCopy, ASTFunctionCall.class)

src/main/java/org/nest/codegeneration/sympy/SympyScriptGenerator.java

@@ -138,22 +138,15 @@ private static Path generateSympyScript(
       Log.warn("It works only for a single ODE. Only the first equation will be used.");
     }
 
-    final ASTEquation workingVersion = ODETransformer.replaceSumCalls(astOdeDeclaration.getODEs().get(0));
-
-    glex.setGlobalValue("ode", workingVersion);
+    glex.setGlobalValue("ode", astOdeDeclaration.getODEs().get(0));
     glex.setGlobalValue("shapes", astOdeDeclaration.getShapes());
     glex.setGlobalValue("predefinedVariables", PredefinedVariables.gerVariables());
 
     setup.setGlex(glex);
     setup.setCommentStart(Optional.of("#"));
     setup.setCommentEnd(Optional.empty());
 
-    final GeneratorEngine generator = new GeneratorEngine(setup);
-
-    final Path solverSubPath = Paths.get( neuron.getName() + "Solver.py");
-
     final Set<VariableSymbol> variables = new HashSet<>(getVariableSymbols(astOdeDeclaration));
-
     final List<VariableSymbol> aliases = ASTUtils.getAliasSymbols(astOdeDeclaration);
 
     List<VariableSymbol> symbolsInAliasDeclaration = aliases
@@ -174,7 +167,10 @@ private static Path generateSympyScript(
 
     final ExpressionsPrettyPrinter expressionsPrinter  = new ExpressionsPrettyPrinter();
     glex.setGlobalValue("printer", expressionsPrinter);
+    glex.setGlobalValue("odeTransformer", new ODETransformer());
 
+    final GeneratorEngine generator = new GeneratorEngine(setup);
+    final Path solverSubPath = Paths.get( neuron.getName() + "Solver.py");
     generator.generate(templateName, solverSubPath, astOdeDeclaration);
 
     return Paths.get(setup.getOutputDirectory().getPath(), solverSubPath.toString());

src/main/resources/org/nest/nestml/neuron/NeuronClass.ftl

@@ -215,7 +215,6 @@ void
 ${simpleNeuronName}::calibrate()
 {
   B_.logger_.init();
-  B_.receptor_types_.resize(1, 0);
 
   <#list body.getInternalNonAliasSymbols() as variable>
     ${tc.includeArgs("org.nest.nestml.function.Calibrate", [variable])}

src/main/resources/org/nest/nestml/neuron/NeuronHeader.ftl

@@ -146,11 +146,11 @@ public:
   <#list body.getInternalSymbols() as internal>
   ${tc.includeArgs("org.nest.nestml.function.MemberVariableGetterSetter", [internal])}
   </#list>
-
+  <#-- DO I NEED THEM?
   <#list body.getODEAliases() as odeAlias>
     ${tc.includeArgs("org.nest.nestml.function.MemberVariableGetterSetter", [odeAlias])}
   </#list>
-
+  -->
   <#list body.getInputBuffers() as buffer>
   ${bufferHelper.printBufferGetter(buffer, false)};
   </#list>

src/main/resources/org/nest/sympy/DeltaShapeSolver.ftl

@@ -8,11 +8,11 @@ __a__, __h__, delta = symbols('__a__ __h__ delta')
 
 # Handle aliases
 <#list aliases as alias>
-${alias.getName()} = ${printer.print(alias.getDeclaringExpression().get())}
+${alias.getName()} = ${printer.print(odeTransformer.replaceSumCalls(alias.getDeclaringExpression().get()))}
 </#list>
 
 # Shapes must be symbolic for the differetiation step
-rhsTmp = ${printer.print(ode.getRhs())}
+rhsTmp = ${printer.print(odeTransformer.replaceSumCalls(ode.getRhs()))}
 constantInputs = simplify(1/diff(rhsTmp, ${shapes[0].getLhs()}) * (rhsTmp - diff(rhsTmp, ${ode.getLhs().getSimpleName()})*${ode.getLhs().getSimpleName()}) - (
 <#assign operator = "">
 <#compress> <#list shapes as eq>
@@ -23,9 +23,9 @@ ${operator} ${eq.getLhs()}
 
 # TODO take the comment for the source model
 <#list shapes as eq>
-${eq.getLhs()} = ${printer.print(eq.getRhs())}
+${eq.getLhs()} = ${printer.print(odeTransformer.replaceSumCalls(eq.getRhs()))}
 </#list>
-rhs = ${printer.print(ode.getRhs())}
+rhs = ${printer.print(odeTransformer.replaceSumCalls(ode.getRhs()))}
 dev${ode.getLhs().getSimpleName()} = diff(rhs, ${ode.getLhs().getSimpleName()})
 dev_t_dev${ode.getLhs().getSimpleName()} = diff(dev${ode.getLhs().getSimpleName()}, t)
 
@@ -39,6 +39,6 @@ if dev_t_dev${ode.getLhs().getSimpleName()} == 0:
     c1 = diff(rhs, ${ode.getLhs().getSimpleName()})
     # The symbol must be declared again. Otherwise, the right hand side will be used for the derivative
     ${shapes[0].getLhs()} = symbols("${shapes[0].getLhs()}")
-    c2 = diff( ${printer.print(ode.getRhs())} , ${shapes[0].getLhs()})
+    c2 = diff( ${printer.print(odeTransformer.replaceSumCalls(ode.getRhs()))} , ${shapes[0].getLhs()})
     f = open('P30.tmp', 'w')
     f.write("P30 real = " + str(simplify(c2 / c1 * (exp(__h__ * c1) - 1))) + "# P00 expression")

src/main/resources/org/nest/sympy/ODESolver.ftl

@@ -8,11 +8,11 @@ __a__, __h__ = symbols('__a__ __h__')
 
 # Handle aliases
 <#list aliases as alias>
-${alias.getName()} = ${printer.print(alias.getDeclaringExpression().get())}
+${alias.getName()} = ${printer.print(odeTransformer.replaceSumCalls(alias.getDeclaringExpression().get()))}
 </#list>
 
 # Shapes must be symbolic for the differetiation step
-rhsTmp = ${printer.print(ode.getRhs())}
+rhsTmp = ${printer.print(odeTransformer.replaceSumCalls(ode.getRhs()))}
 constantInputs = simplify(1/diff(rhsTmp, ${shapes[0].getLhs()}) * (rhsTmp - diff(rhsTmp, ${ode.getLhs().getSimpleName()})*${ode.getLhs().getSimpleName()}) - (
 <#assign operator = "">
 <#compress> <#list shapes as eq>
@@ -23,9 +23,9 @@ ${operator} ${eq.getLhs()}
 
 # TODO take the comment for the source model
 <#list shapes as eq>
-${eq.getLhs()} = ${printer.print(eq.getRhs())}
+${eq.getLhs()} = ${printer.print(odeTransformer.replaceSumCalls(eq.getRhs()))}
 </#list>
-rhs = ${printer.print(ode.getRhs())}
+rhs = ${printer.print(odeTransformer.replaceSumCalls(ode.getRhs()))}
 dev${ode.getLhs().getSimpleName()} = diff(rhs, ${ode.getLhs().getSimpleName()})
 dev_t_dev${ode.getLhs().getSimpleName()} = diff(dev${ode.getLhs().getSimpleName()}, t)
 
@@ -99,7 +99,7 @@ if dev_t_dev${ode.getLhs().getSimpleName()} == 0:
         c1 = diff(rhs, ${ode.getLhs().getSimpleName()})
         # The symbol must be declared again. Otherwise, the right hand side will be used for the derivative
         ${shapes[0].getLhs()} = symbols("${shapes[0].getLhs()}")
-        c2 = diff( ${printer.print(ode.getRhs())} , ${shapes[0].getLhs()})
+        c2 = diff( ${printer.print(odeTransformer.replaceSumCalls(ode.getRhs()))} , ${shapes[0].getLhs()})
 
         # define matrices depending on order
         # for order 1 and 2 A is lower triangular matrix

src/test/java/org/nest/integration/NESTMLFrontendIntegrationTest.java

@@ -49,7 +49,7 @@ public void testInfrastructure() {
   @Test
   public void testManually() {
     final String[] args = new String[] {
-        "models/aeif_cond_alpha_implicit.nestml",
+        "models/iaf_psc_exp.nestml",
         "--target", outputPath.toString()};
 
     new NESTMLFrontend().start(args);

src/test/java/org/nest/integration/ODEProcessorTest.java

@@ -67,7 +67,6 @@ public void testDeltaModel() throws Exception {
   private Scope processModel(final String pathToModel) {
     final ASTNESTMLCompilationUnit modelRoot = parseNESTMLModel(pathToModel);
     scopeCreator.runSymbolTableCreator(modelRoot);
-    final String modelFolder = modelRoot.getFullName();
     final Path outputBase = Paths.get(OUTPUT_FOLDER.toString(), Names.getPathFromQualifiedName(pathToModel.toString()));
     FilesHelper.deleteFilesInFolder(outputBase);