Hodgkin-Huxley model in C++ I'm fairly new to C++ and I want to simulate the Hodgkin-Huxley neuron model with it. I have used a MATLAB implementation before and I hope that the C++ code will be faster.  

Both seem to take the same time when I time them. For the C++ implementation I used the Linux time command and in MATLAB tic/toc. I let both simulations run for 1000 times with 100000 time steps each. Both take around 13-15s to complete.  

I hope anybody can tell me if I have some mayor issues in my C++ implementation or if the difference between C++/MATLAB for such a task is just small.   

I compiled the C++ code with g++ -std=gnu++11 -O3 -o hh hhModel.cpp

#include &lt;iostream&gt;
#include &lt;cstdlib&gt;
#include &lt;stdio.h&gt;
#include &lt;math.h&gt;
#include &lt;array&gt;
using namespace std;

int main()
{
    double dt = 0.01;
    int numSamples = 100000;

    // create an rectangular external stimulus
    double I[numSamples];
    for(int i = 0; i &lt;= numSamples; i++)
        I[i] = 0;
    for(int i = 25000; i &lt;= 75000; i++)
            I[i] = 0.01;

    // init constants
    double Vinit = -65;
    double Vref = -65;

    double Smemb = 4000;  // [um^2] surface area of the membrane
    double Cmemb = 1.0;   // [uF/cm^2] membrane capacitance density
    double Cm = Cmemb * Smemb * 1e-8;   // [uF] membrane capacitance

    double GNa = 120;
    double GK = 36;
    double GL = 0.3;
    double ENa = 125;
    double EK = -55;
    double EL = -25;

    double v[numSamples]; // [mV] membrane potential vector

    double m, h, n, aM, bM, aH, bH, aN, bN, dv_dt, dn_dt, dm_dt, dh_dt, IL, IK, INa;

    double gNa = GNa * Smemb * 1e-8;// Na conductance [mS]
    double gK = GK * Smemb * 1e-8;  // K conductance [mS]
    double gL = GL * Smemb * 1e-8;  // leak conductance [mS]

    // initial values
    v[0] = Vinit;           // initial membrane potential
    m = 0;
    h = 0;
    n = 0;

    // To compare the code to Matlab I run the numerical integration 1000 times
    for (int nn = 0; nn &lt;= 1000; nn++){

    // Numerical integration
    for (int j = 0; j &lt;= numSamples; j++){
        // ionic currents
        INa = gNa * m*m*m * h * (ENa - v[j]);
        IK = gK * n*n*n*n * (EK - v[j]);
        IL = gL * (EL - v[j]);

        aM = 0.1 * (v[j] - Vref -25) / ( 1 - exp(-(v[j]-Vref-25)/10));
        bM = 4 * exp(-(v[j]-Vref)/18);

        aH = 0.07 * exp(-(v[j]-Vref)/20);
        bH = 1 / ( 1 + exp( -(v[j]-Vref-30)/10 ) );

        aN = 0.01 * (v[j]-Vref-10) / ( 1 - exp(-(v[j]-Vref-10)/10) );
        bN = 0.125 * exp(-(v[j]-Vref)/80);

        // derivatives
        dv_dt = ( INa + IK + IL + I[j] ) / Cm;
        dm_dt = (1.0-m) * aM - m * bM;
        dh_dt = (1.0-h) * aH - h * bH;
        dn_dt = (1.0-n) * aN - n * bN;

        // calculate next step
        v[j+1] = v[j] + dv_dt * dt;
        m = m + dm_dt * dt;
        h = h + dh_dt * dt;
        n = n + dn_dt * dt;
    }
    }
}
Note: benchmarking C++ code over a fixed, compiler provided sample can lead to the optimizer removing all relevant code and/or replacing it by a precalculated result. The only reason this didn't happen in your test was that GCC gave up because of the complexity and not enough hints.


Implementation


I and v are declared with auto storage class. This causes them to be placed on the stack on most platforms, for which they might be too large (2 * 100000 * 8 Bytes = 1.6 MB, for comparison: standard max stack size on windows is 1 MiB). This might crash the program or cause any number of undefined behaviors.
Also, double I[numSamples]; and double v[numSamples]; would require numSamples to be a compile time constant (e.g. const or constexpr) by the C++ standard (GCC seems not to mind, though).
Checking array/loop bounds: All loop bounds are off by one (&lt;= instead of &lt;), which in turn causes array accesses to be out of bounds by 1 or even 2. This is undefined behavior, and can lead to crashes! (Also, it means you've run 1001 loops of the C++ code, instead of the 1000 mentioned, slightly skewing the comparison.)
A lot of variables could be made const and/or constexpr. Doing so might help the compiler optimize your code, as it gets more information to work with/doesn't have to deduce the const-ness by itself (which can be costly).
Unnecessary header #includes: contents from &lt;iostream&gt;, &lt;stdio.h&gt; and &lt;array&gt; are never used!
Prefer to include the adapted C++ headers over plain C headers (&lt;cmath&gt; and &lt;cstdio&gt; instead of &lt;math.h&gt; and &lt;stdio.h&gt;). 
using namespace std; is considered bad practice and should be avoided.
Try to refactor code into smaller, independent units. This helps humans for readability and compilers for deducing useful traits for optimization.
Along the same line: Try to restrict variables to the smallest scope needed. Again, this helps human readers as much as compilers, as now we only have to consider them where they are reachable.
Your code still uses a lot of magic numbers. Try to extract them into appropriately named variables, so other readers of your code can reason about them easier.


Optimization

Your current code is not structured well enough for the optimizer. Why? If it could reason about it perfectly, it would remove all the code (as there are no observable side effects), running in near 0ms!

Declaring all possible variables const and moving them to the smallest scope possible helps quite a bit. Even better: Refactor the code into functions (and maybe classes) with proper const annotations!

Just this alone allowed MSVC and GCC to recognize that the whole calculations aren't needed (as they don't cause side effects by themselves), so they removed them and run-time went down to ~24 microseconds on my machine.

To get relevant results, we now need to fool the compiler a bit: Adding a little side effect (num_calls in the code below) and taking the number of loops to run as a command line argument were enough for MSVC and GCC (and clang, though not ICC). Now runtime (with all calculations!) went down from ~14 seconds (original) to ~12.5 seconds (again on my machine).


  Note: The code below has some changes that, while possibly benefiting benchmark performance, are differences to the original:
  
  
  Results of all the simulation steps aren't stored. Instead, only the latest result is kept. (The stores would probably be optimized out anyways, as they wouldn't be read from.)
  The input impulse is generated on the fly (instead of being read from an array). This causes the optimizer to generate better code for this case (no memory reads, just a branch that the hardware can predict well). I don't think the impact is large (as contiguous array accesses would be prefetched anyways), but still to be noted.
  The number of benchmark loops to run is taken from the command line arguments (argument #1). This is done to fool the optimizer into not removing/precalculating the actual calls (as it can't know how often it will run).
  To cause a visible side effect, I needed to report the number of simulation steps made. Since I had to print them anyways, I added a small snippet to measure the benchmark time and print it as well.
  I renamed a lot of variables. This was more for my understanding (and readability). For some variables, I couldn't deduce their intents, so their names remained unchanged. I also didn't move the magic numbers into variables, as I don't know what they are referring to. (I know scientific code usually tries to use formula symbols as names, but those were too similar for me to reason about, especially since I'm not too versed into that specific field.)
  


#include &lt;iostream&gt;
#include &lt;cmath&gt;
#include &lt;chrono&gt;

uint32_t num_calls = 0; // to fool optimizer

namespace constants
{
    constexpr const auto time_step = 0.01;
    constexpr const auto num_samples = 100000;

    constexpr const auto initial_voltage = -65.0;
    constexpr const auto reference_voltage = -65.0;

    constexpr const auto membrane_surface = 4000.0;  // [um^2] surface area of the membrane
    constexpr const auto membrane_capacitance_density = 1.0;   // [uF/cm^2] membrane capacitance density
    constexpr const auto membrane_capacitance = membrane_capacitance_density * membrane_surface * 1e-8;   // [uF] membrane capacitance

    constexpr const auto GNa = 120.0;
    constexpr const auto GK = 36.0;
    constexpr const auto GL = 0.3;
    constexpr const auto ENa = 125.0;
    constexpr const auto EK = -55.0;
    constexpr const auto EL = -25.0;

    constexpr const auto sodium_conductance = GNa * membrane_surface * 1e-8;// Na conductance [mS]
    constexpr const auto kalium_conductance = GK * membrane_surface * 1e-8;  // K conductance [mS]
    constexpr const auto leak_conductance = GL * membrane_surface * 1e-8;  // leak conductance [mS]
}

class simulation
{
    double membrane_potential;
    double m;
    double h;
    double n;

    double step_m(const double delta_time) const
    {
        const auto voltage_difference = membrane_potential - constants::reference_voltage;

        const auto aM = 0.1 * (voltage_difference - 25.0) / (1.0 - exp(-(voltage_difference - 25.0) / 10.0));
        const auto bM = 4.0 * exp(-voltage_difference / 18.0);

        auto const derivative_m = (1.0 - m) * aM - m * bM;

        return (m + derivative_m * delta_time);
    }

    double step_h(const double delta_time) const
    {
        const auto voltage_difference = membrane_potential - constants::reference_voltage;

        const auto aH = 0.07 * exp(-voltage_difference / 20.0);
        const auto bH = 1.0 / (1.0 + exp(-(voltage_difference - 30.0) / 10.0));

        auto const derivative_h = (1.0 - h) * aH - h * bH;

        return h + derivative_h * delta_time;
    }

    double step_n(const double delta_time) const
    {
        const auto voltage_difference = membrane_potential - constants::reference_voltage;

        const auto aN = 0.01 * (voltage_difference - 10.0) / (1.0 - exp(-(voltage_difference - 10.0) / 10.0));
        const auto bN = 0.125 * exp(-voltage_difference / 80.0);

        auto const derivative_n = (1.0 - n) * aN - n * bN;

        return n + derivative_n * delta_time;
    }

    double step_membrane_potential(const double time, const double stimulus) const
    {
        const auto sodium_current = constants::sodium_conductance * m*m*m * h * (constants::ENa - membrane_potential);
        const auto kalium_current = constants::kalium_conductance * n*n*n*n * (constants::EK - membrane_potential);
        const auto leak_current = constants::leak_conductance * (constants::EL - membrane_potential);

        auto const derivative_potential = (sodium_current + kalium_current + leak_current + stimulus) / constants::membrane_capacitance;

        return membrane_potential + derivative_potential * time;
    }
public:
    simulation(const double initial_potential, const double initial_m, const double initial_h, const double initial_n) : membrane_potential(initial_potential), m(initial_m), h(initial_h), n(initial_n) {}

    void step(const double delta_time, const double current_stimulus)
    {
        membrane_potential = step_membrane_potential(delta_time, current_stimulus);
        m = step_m(delta_time);
        h = step_h(delta_time);
        n = step_n(delta_time);

        ++num_calls; // to fool optimizer
    }

    double latest_potential() const
    {
        return membrane_potential;
    }
};

void run_simulation()
{
    auto sim = simulation{ constants::initial_voltage, 0, 0, 0 };

    for (auto j = 0; j &lt; constants::num_samples; j++) {
        auto stimulus = 0.0;
        if (j &gt;= 25000 &amp;&amp; j &lt;= 75000) stimulus = 0.01;

        sim.step(constants::time_step, stimulus);
    }
}

int main(int argc, char *argv[])
{
    using seconds = std::chrono::duration&lt;double, std::ratio&lt;1, 1&gt;&gt;;

    if(argc &lt; 2)
    {
        std::cerr &lt;&lt; "requires number of benchmarking loops as command line argument #1";
        return 1;
    }

    const auto num_loops = atoi(argv[1]);

    const auto start_time = std::chrono::high_resolution_clock::now();

    for (auto loop = 0; loop &lt; num_loops; ++loop) {
        run_simulation();
    }

    const auto end_time = std::chrono::high_resolution_clock::now();
    const auto time_taken = std::chrono::duration_cast&lt;seconds&gt;(end_time - start_time);

    std::cout &lt;&lt; "time taken: " &lt;&lt; time_taken.count() &lt;&lt; "s\n";
    std::cout &lt;&lt; "number of loops: " &lt;&lt; num_loops &lt;&lt; "\n";
    std::cout &lt;&lt; "number of simulation steps: " &lt;&lt; (num_calls / num_loops) &lt;&lt; "\n"; // to fool optimizer
}


Further optimization possibilities?


All compilers produce reasonable ILP (instruction level parallelism) assembly (if still fooled). There doesn't seem to be much to improve upon.
Multithreading also doesn't seem like it would be helping for this problem, as the simulation steps are dependent on previous results.
The bottleneck seems to be the exp function call, so a faster implementation might help (that's why I tried ICC, but that one is too smart and removes the calculations as a whole -.-).