The C++ code below simulates the timecourse of the membrane potential (V) of a population of 128 leaky integrate-and-fire neurons using the Euler method for numerical integration.

It produces data for three data points in Figure 1 of this paper. (Well, it produces 400 ms worth of data for three data points.)

I have an equivalent script written in Python (which I'm happy to share, if anyone wants to help optimize that). Expecting a major speed boost with C++ (e.g., 10-100X), I was surprised to see that this script took 20.7 s to run on my laptop, as compared to 77.5 s with Python (a smaller than 4X boost).

But, since I'm a C++ newbie, I'm hoping there's work I can do to optimize this. (I'd also love to get style critiques, because I'm sure this is stylistically a mess.)

This script saves three text files to a directory on my laptop. In order to run this, you will need to change the directory specification to a directory that exists on your computer.

#include <math.h>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <iterator>
#include <Eigen/Dense>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
using namespace std;
using Eigen::MatrixXd;
using Eigen::ArrayXd;

bool save_mat(MatrixXd pdata, const string& file_path)
  ofstream os(file_path.c_str());
  if (!os.is_open())
      cout << "Failure!" << endl;
      return false;
  int n_rows = pdata.rows();
  int n_cols = pdata.cols();
  for (int i = 0; i < n_rows; i++)
      for (int j = 0; j < n_cols; j++)
          os << pdata(i, j);
          if (j + 1 == n_cols)
              os << endl;
              os << ",";
  return true;

template<typename T>
vector<T> arange(T start, T stop, T step = 1)
  vector<T> values;
  for (T value = start; value < stop; value += step)
  return values;

double f(double t, double tau_1, double tau_2);

ArrayXd set_initial_V(double tau, double g_L, double I_0, double theta,
                      double V_L, int N, double c);

int main()

  double tau_1, tau_2, I_syn_bar, dt, tau, g_L, I_0, theta, V_L, c, t, C;
  int N, n_t;
  string dan_dir;

  tau_1 = 3.0;  // ms.
  tau_2 = 1.0;  // ms.
  I_syn_bar = 0.5; // uA / cm^2.
  N = 128;
  double dt_array[3] = {0.25, 0.1, 0.01};  // ms.
  tau = 10;  // ms.
  g_L = 0.1;  // mS / cm^2.
  I_0 = 2.3;  // uA / cm^2.
  theta = -40;  // mV.
  V_L = -60;  // mV.
  c = 0.5;
  C = 1;  // uF / cm^2.

  struct stat statbuf;
  if (stat("/home/daniel", &statbuf) == 0 && S_ISDIR(statbuf.st_mode))
      dan_dir = "/home/daniel/Science";
  else if (stat("/home/dan", &statbuf) == 0 && S_ISDIR(statbuf.st_mode))
      dan_dir = "/home/dan/Science";
  string save_file = ("/dopa_net/results/hansel/test/test_hansel%.2g.txt");
  save_file = dan_dir + save_file;
  char complete_save_file[save_file.length()];

  int n_dt = (sizeof(dt_array) / sizeof(*dt_array));

  for (int k = 0; k < n_dt; k++)
      dt = dt_array[k];
      cout << dt << endl;
      auto t_vector = arange<double>(0, 400, dt);
      n_t = t_vector.size();
      MatrixXd V(N, n_t);
      V.col(0) = set_initial_V(tau, g_L, I_0, theta, V_L, N, c);
      double I_syn = 0;  // uA / cm^2.
      ArrayXd t_spike_array = ArrayXd::Zero(N);
      for (int i = 1; i < n_t; i++)
          auto prev_V = V.col(i - 1).array();
          auto current_V = prev_V + dt * (-g_L * (prev_V - V_L) + I_syn +
                                          I_0) / C;
          V.col(i) = current_V;
          t = t_vector[i];
          I_syn = 0;
          for (int j = 0; j < N; j++)
              if (current_V(j) > theta)
                  t_spike_array(j) = t;
                  V(j, i) = V_L;
              I_syn += I_syn_bar / N * f(t - t_spike_array(j), tau_1, tau_2);
      sprintf(complete_save_file, save_file.c_str(), dt);
      save_mat(V, complete_save_file);
  return 0;


double f(double t, double tau_1, double tau_2)
  return tau_2 / (tau_1 - tau_2) * (exp(-t / tau_1) - exp(-t / tau_2));

ArrayXd set_initial_V(double tau, double g_L, double I_0, double theta,
                      double V_L, int N, double c)
  double T;
  int i;
  T = -tau * log(1 - g_L / I_0 * (theta - V_L));
  ArrayXd V(N);
  for (i = 0; i < N; i++)
      V(i) = V_L + I_0 / g_L * (1 - exp(-c * (i - 1) / N * T / tau));
  return V;
  • \$\begingroup\$ i just got rid of the C++11 dependency by removing the use of auto. this may or may not have given me a significant speed boost. on a different computer from the one i used for the timing info in the question, i'm getting a run time of 8.8 s with C++. \$\endgroup\$
    – abcd
    Commented Oct 15, 2015 at 7:34
  • \$\begingroup\$ Removing use of auto will have no effect on performance. You could always profile the code to find out exactly what takes up the time. Also, investigate your compilers optimization settings. \$\endgroup\$
    – user673679
    Commented Oct 15, 2015 at 11:43

1 Answer 1



Compiler flags

Your main performance gain will simple be using the -O3 flag for compilation. You may very well achieve even better performance with more sophisticated flags depending on your compiler and system. On my system it takes ~4s.


Using a the simple linux profiling tool perf, shows that > 60 % of the remaining time is actually spent in vfprintf. A specific performance analysis tool recommendation depends on your OS.

Disabling the output, the remaining runtime is ~0.2s.

So you could now go ahead and try to optimize the output formatting, but is it really worth it for your use case?


Don't use sprintf

sprintf is very dangerous. Don't use it, not in C, nor in C++. In C++ concatenate your std::strings with + or use std::stringstream. If you really want format like syntax use boost::format.

Split up your main function

Don't handle file name computations and math stuff in the same function. It's very distracting.

Dont use using namespace std;

Yes, it requires more typing, but its bad practice. For an explanation take a look at the answers to this question.

Use (const) references when appropriate

Disclaimer: I don't know Eigen, but I assume MatrixXd is an expensive-to-copy object.

You should use references to hand expensive-to-copy objects to functions. If they are input only, use const references - if they are in/out non-const references. So for save_mat use a const MatrixXd& just like you do for the file name string.

One could argue that auto current_V = prev_V ... is also dangerously expensive because it requries assignment of expensive-to-copy objects. But in this case it is an rvalue so move assignment can be used. You could also use const auto& current V = ... instead which would express the intent of not modifying it later and perform slightly better (the latter being probably irrelevant)

Note The rest is mainly C++11 biased.


In C++11, avoid raw loops. Use range-based loops instead.

Especially with the vector rcreated by arrange - It's neither efficient nor very clear to store the values you are going to process in the t_vector that you can easily generate on the fly instead of wasting memory on that. Usually I would recommend looking at cppitertools or range-v3, which allows for a more efficient and cleaner abstraction of the t values.

Also use range-based loops for the outer loop through dt_array. You can also replace the innermost loop with with cppitertools/range-v3.

Use const or constexpr....

... for local variables that don't change. Further define one variable per line and set it's value right at the definition (even if it is not const).

Use std::array instead of array[]

Simply use a constexpr n_dt and stay away for those sizeof hacks to figure out the number of elements you just defined a few lines earlier.

  • \$\begingroup\$ Usually compute clusters provide multiple compilers in multiple versions e.g. via gnu modules. Please consult the documentation or administration about what compilers are actually available and if new versions can be installed. It is absolutely worth it to use C++11/14, especially as a beginner. Modern C++ is far superior in terms of clarity and ease of safe and performant use. \$\endgroup\$
    – Zulan
    Commented Oct 16, 2015 at 7:53
  • \$\begingroup\$ Does your ~5 hour problem really have the same output / computational ratio? \$\endgroup\$
    – Zulan
    Commented Oct 16, 2015 at 15:02
  • \$\begingroup\$ What I'm saying: The code code generates ~80 MB in 4 seconds (IIRC). Considering the same output / computation ratio, this would be 360 GB in 5 hours - right? What kind of post-processing do you do with those 360 GB - and how does the post-processing time relate to the 5 hours of data generation. \$\endgroup\$
    – Zulan
    Commented Oct 16, 2015 at 15:11
  • \$\begingroup\$ my OS (CentOS6/RHEL6) doesn't have perf. is there another specific tool you'd recommend i use? \$\endgroup\$
    – abcd
    Commented Oct 16, 2015 at 17:22
  • \$\begingroup\$ new question posted here. \$\endgroup\$
    – abcd
    Commented Oct 16, 2015 at 19:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.