Simple neural-network simulation in C++

Question

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;
    }
  os.precision(11);
  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;
            }
          else
            {
              os << ",";
            }
        }
    }
  os.close();
  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)
    values.push_back(value);
  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;
}

Answer 1

Performance

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.

Output

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?

Style

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.

Loops

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.