Merge pull request #140 from DimitriPlotnikov/master

Make gif models more consistent. Fix the ODE analyser: use GSL sovler in case of many ODEs in 'normal' form.

src/main/java/org/nest/codegeneration/NESTCodeGenerator.java

@@ -66,7 +66,7 @@ public void analyseAndGenerate(
 
     ASTNESTMLCompilationUnit workingVersion = root;
     for (int i = 0; i < root.getNeurons().size(); ++i) {
-      final ASTNeuron solvedNeuron = computeSolutionForODE(root.getNeurons().get(i), outputBase);
+      final ASTNeuron solvedNeuron = solveODESInNeuron(root.getNeurons().get(i), outputBase);
       root.getNeurons().set(i, solvedNeuron);
     }
 
@@ -81,7 +81,7 @@ public void analyseAndGenerate(
     info(msg, LOG_NAME);
   }
 
-  private ASTNeuron computeSolutionForODE(
+  private ASTNeuron solveODESInNeuron(
       final ASTNeuron astNeuron,
       final Path outputBase) {
     final ASTBody astBody = astNeuron.getBody();
@@ -91,8 +91,11 @@ private ASTNeuron computeSolutionForODE(
         info("The model will be solved numerically with a GSL solver.", LOG_NAME);
         return astNeuron;
       }
+      else {
+        info("The model will be analysed.", LOG_NAME);
+        return odeProcessor.solveODE(astNeuron, outputBase);
+      }
 
-      return odeProcessor.solveODE(astNeuron, outputBase);
     }
     else {
       return astNeuron;
@@ -294,12 +297,11 @@ private void setNeuronGenerationParameter(
   }
 
   private void defineSolverType(final GlobalExtensionManagement glex, final ASTNeuron neuron) {
-    Log.warn("This logic doesn't work at the moement. The sympy script must be used to consult which ODE how should be solved.");
     final ASTBody astBody = neuron.getBody();
     glex.setGlobalValue("useGSL", false);
     if (astBody.getEquations().isPresent()) {
-      if (astBody.getEquations().get().getODEs().size() > 1) {
-        glex.setGlobalValue("useGSL", false);
+      if (astBody.getEquations().get().getEqs().size() == 0) {
+        glex.setGlobalValue("useGSL", true);
         glex.setGlobalValue("ODEs", astBody.getEquations().get().getODEs());
       }
 

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

@@ -37,7 +37,6 @@
 import static java.util.Optional.of;
 import static java.util.stream.Collectors.toList;
 import static org.nest.symboltable.predefined.PredefinedFunctions.I_SUM;
-import static org.nest.symboltable.predefined.PredefinedFunctions.SUM;
 import static org.nest.utils.ASTNodes.getVariableSymbols;
 
 /**
@@ -98,7 +97,7 @@ private static Path generateSolverScript(
       Log.warn("It works only for a single ODE. Only the first equation will be used.");
     }
     ASTODE workingVersion = replace_I_sum(astOdeDeclaration.getODEs().get(0));
-    workingVersion = replace_sum(workingVersion);
+
     glex.setGlobalValue("ode", workingVersion);
     glex.setGlobalValue("EQs", astOdeDeclaration.getEqs());
     glex.setGlobalValue("predefinedVariables", PredefinedVariables.gerVariables());
@@ -157,21 +156,6 @@ static ASTODE replace_I_sum(final ASTODE astOde) {
     return astOde;
   }
 
-  private static ASTODE replace_sum(final ASTODE astOde) {
-    final List<ASTFunctionCall> functions = ASTNodes.getAll(astOde, ASTFunctionCall.class)
-        .stream()
-        .filter(astFunctionCall -> astFunctionCall.getCalleeName().equals(SUM))
-        .collect(toList());
-
-    functions.stream().forEach(node -> {
-      final Optional<ASTNode> parent = ASTNodes.getParent(node, astOde);
-      checkState(parent.isPresent());
-      final ASTExpr expr = (ASTExpr) parent.get();
-      expr.setFunctionCall(null);
-      expr.setVariable(node.getArgs().get(0).getVariable().get());
-    });
-    return astOde;
-  }
 
   private static GlobalExtensionManagement createGLEXConfiguration() {
     return new GlobalExtensionManagement();

src/main/java/org/nest/symboltable/predefined/PredefinedFunctions.java

@@ -36,7 +36,6 @@ public class PredefinedFunctions {
   private static final String RANDOM_INT = "randomInt";
   private static final String EXPM1 = "expm1";
   private static final String DELTA = "delta";
-  public static final String SUM = "sum";
   private static final String MAX = "max";
   public static final String INTEGRATE = "integrate";
   public static final String I_SUM = "I_sum";
@@ -108,11 +107,6 @@ public class PredefinedFunctions {
     delta.setReturnType(getType("real"));
     name2FunctionSymbol.put(DELTA, delta);
 
-    final MethodSymbol sum = createFunctionSymbol(SUM);
-    sum.addParameterType(getType("pA"));
-    sum.setReturnType(getType("pA"));
-    name2FunctionSymbol.put(SUM, sum);
-
     final MethodSymbol max = createFunctionSymbol(MAX);
     max.addParameterType(getType("mV"));
     max.addParameterType(getType("mV"));

src/test/resources/codegeneration/gif/glif.nestml

@@ -3,16 +3,16 @@ neuron glif:
   state:
     V mV
     I_a pA[n]
-    I_shape_in, I_shape_ex real
+    I_shape_in, I_shape_ex, I_shape_a real
     Theta   mV
   end
 
   equations:
-      I_shape_in = (e/tau_syn_in) * t * exp(-1/tau_syn_in*t)
-      I_shape_ex = (e/tau_syn_ex) * t * exp(-1/tau_syn_ex*t)
-      I_a' = -(1/tau_a)*I_a
+      I_shape_in = exp(-1/tau_syn_in*t)
+      I_shape_ex = exp(-1/tau_syn_ex*t)
+      I_shape_a =  exp(-1/tau_syn_a*t)
 
-      V' = -1/Tau * (V-E_L) + 1/C_m * (sum(I_a)+I_sum(I_shape_in, in_spikes) + I_sum(I_shape_ex, ex_spikes) + I_e + currents)
+      V' = -1/Tau * (V-E_L) + 1/C_m * (I_shape_a+I_sum(I_shape_in, in_spikes) + I_sum(I_shape_ex, ex_spikes) + I_e + currents)
 
       Theta' = -b*(Theta-Theta_inf) + a*(V-E_L)
   end
@@ -27,7 +27,7 @@ neuron glif:
     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.
-    tau_a ms[n] = 2     # Time constant of synaptic current.
+    tau_syn_a ms[n] = 2     # Time constant of synaptic current.
     t_ref   ms = 2     # Refractory period.
     E_L     mV = -70  # Resting potential.
     V_reset mV = 0 [[V_reset < Theta_reset]]

src/test/resources/codegeneration/gif/glif_2.nestml

@@ -9,8 +9,8 @@ neuron glif_2:
   end
 
   equations:
-      I_shape_in = (e/tau_syn_in) * t * exp(-1/tau_syn_in*t)
-      I_shape_ex = (e/tau_syn_ex) * t * exp(-1/tau_syn_ex*t)
+      I_shape_in = exp(-1/tau_syn_in*t)
+      I_shape_ex = exp(-1/tau_syn_ex*t)
 
       Theta' = -b*(Theta-Theta_inf) + a*(V-E_L)
       I_a_1' = -(1/tau_a_1)*I_a_1

src/test/resources/codegeneration/gif/glif_extended.nestml

@@ -8,12 +8,12 @@ neuron glif_extended:
   end
 
   equations:
-      I_shape_in = (e/tau_syn_in) * t * exp(-1/tau_syn_in*t)
-      I_shape_ex = (e/tau_syn_ex) * t * exp(-1/tau_syn_ex*t)
+      I_shape_in = exp(-1/tau_syn_in*t)
+      I_shape_ex = exp(-1/tau_syn_ex*t)
 
       Theta' = -b*(Theta-Theta_inf) + a*(V-E_L)
       I_a' = Beta*(V-E_L) - (1/tau_a)*I_a
-      V' = -1/Tau * (V-E_L) + 1/C_m * (I_sum(I_shape_in, in_spikes) + I_sum(I_shape_ex, ex_spikes) + I_e + currents + sum(I_a))
+      V' = -1/Tau * (V-E_L) + 1/C_m * (I_sum(I_shape_in, in_spikes) + I_sum(I_shape_ex, ex_spikes) + I_e + currents + I_a)
   end
 
   parameter:

src/test/resources/codegeneration/iaf_cond_alpha_implicit.nestml

@@ -17,23 +17,23 @@ neuron iaf_cond_alpha_implicit:
   end
 
   parameter:
-   V_th mV = -55.0      # Threshold Potential in mV
-   V_reset mV = -60.0   # Reset Potential in mV
-   t_ref ms = 2.0       # Refractory period in ms
-   g_L ms = 16.6667     # Leak Conductance in nS
-   C_m pF = 250.0       # Membrane Capacitance in pF
-   alias Tau ms = (1 / g_L) * C_m
-   V_reversalE mV = 0          # Excitatory reversal Potential in mV
-   V_reversalI mV = -85.0      # Inhibitory reversal Potential in mV
-   E_L mV = -70.0       # Leak reversal Potential (aka resting potential) in mV
-   tau_synE ms = 0.2    # Synaptic Time Constant Excitatory Synapse in ms
-   tau_synI ms = 2.0    # Synaptic Time Constant for Inhibitory Synapse in ms
-   I_e pA = 0           # Constant Current in pA
-
-   # Input current injected by CurrentEvent.
-   # This variable is used to transport the current applied into the
-   # _dynamics function computing the derivative of the state vector.
-   I_stim pA = 0
+    V_th mV = -55.0      # Threshold Potential in mV
+    V_reset mV = -60.0   # Reset Potential in mV
+    t_ref ms = 2.0       # Refractory period in ms
+    g_L ms = 16.6667     # Leak Conductance in nS
+    C_m pF = 250.0       # Membrane Capacitance in pF
+    alias Tau ms = (1 / g_L) * C_m
+    V_reversalE mV = 0          # Excitatory reversal Potential in mV
+    V_reversalI mV = -85.0      # Inhibitory reversal Potential in mV
+    E_L mV = -70.0       # Leak reversal Potential (aka resting potential) in mV
+    tau_synE ms = 0.2    # Synaptic Time Constant Excitatory Synapse in ms
+    tau_synI ms = 2.0    # Synaptic Time Constant for Inhibitory Synapse in ms
+    I_e pA = 0           # Constant Current in pA
+
+    # Input current injected by CurrentEvent.
+    # This variable is used to transport the current applied into the
+    # _dynamics function computing the derivative of the state vector.
+    I_stim pA = 0
   end
 
   function set_Tau(v ms):