BlueBrain / nmodl

Code Generation Framework For NEURON MODeling Language
https://bluebrain.github.io/nmodl/
Apache License 2.0
55 stars 15 forks source link

MOD files containing two SOLVEs don't compile. #1501

Open 1uc opened 2 weeks ago

1uc commented 2 weeks ago

For example the following:

$ cat two_blocks.mod
NEURON {
    SUFFIX two_blocks
    RANGE A, B, X, Y
    NONSPECIFIC_CURRENT il
}

STATE {
    X
    Y
    A
    B
}

ASSIGNED {
    il
    i
}

BREAKPOINT {
    SOLVE state1 METHOD sparse
    SOLVE state2 METHOD sparse
    il = i
}

INITIAL {
    A = 1.0
    B = 3.0
    X = 0.0
    Y = 1.0
}

KINETIC state1 {
    ~ X <-> Y (0.4, 0.5)
}

KINETIC state2 {
    ~ A <-> B (0.6, 0.7)
    i = (f_flux-b_flux)
}

generates the following C++ code:

            Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
            double* nmodl_eigen_x = nmodl_eigen_xm.data();
            nmodl_eigen_x[static_cast<int>(0)] = inst->X[id];
            nmodl_eigen_x[static_cast<int>(1)] = inst->Y[id];
            // call newton solver
            functor_two_blocks_0 newton_functor(nt, inst, id, pnodecount, v, indexes, data, thread);
            newton_functor.initialize();
            int newton_iterations = nmodl::newton::newton_solver(nmodl_eigen_xm, newton_functor);
            if (newton_iterations < 0) assert(false && "Newton solver did not converge!");
            inst->X[id] = nmodl_eigen_x[static_cast<int>(0)];
            inst->Y[id] = nmodl_eigen_x[static_cast<int>(1)];
            newton_functor.initialize(); // TODO mimic calling F again.
            newton_functor.finalize();

            Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
            double* nmodl_eigen_x = nmodl_eigen_xm.data();
            nmodl_eigen_x[static_cast<int>(0)] = inst->A[id];
            nmodl_eigen_x[static_cast<int>(1)] = inst->B[id];
            // call newton solver
            functor_two_blocks_1 newton_functor(nt, inst, id, pnodecount, v, indexes, data, thread);
            newton_functor.initialize();
            int newton_iterations = nmodl::newton::newton_solver(nmodl_eigen_xm, newton_functor);
            if (newton_iterations < 0) assert(false && "Newton solver did not converge!");
            inst->A[id] = nmodl_eigen_x[static_cast<int>(0)];
            inst->B[id] = nmodl_eigen_x[static_cast<int>(1)];
            newton_functor.initialize(); // TODO mimic calling F again.
            newton_functor.finalize();

and fails with:

$ nrnivmodl -coreneuron
x86_64/corenrn/mod2c/two_blocks.cpp:786:41: error: redeclaration of ‘Eigen::Matrix<double, 2, 1> nmodl_eigen_xm’
  786 |             Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
      |                                         ^~~~~~~~~~~~~~
x86_64/corenrn/mod2c/two_blocks.cpp:771:41: note: ‘Eigen::Matrix<double, 2, 1> nmodl_eigen_xm’ previously declared here
  771 |             Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
      |                                         ^~~~~~~~~~~~~~

NOCMODL seems to generate code that compiles and looks sane. Therefore, it's not obvious the MOD file is invalid.

1uc commented 2 weeks ago

Summary of discussion with @nrnhines:

  1. Carefully check if NOCMODL supports this and check ModelDB CI for usage examples.
  2. Also consider recursive use, e.g. solving a LINEAR block inside a KINETIC block.
  3. Note that STATE variables are declared globally, i.e. there's no notion of the STATEs of block state1. Some solvers might require knowing the number of states.

Generally, this seems to be a problem that can be revisited when there's an instance of this problem.