Simple neural-network simulation in C++ (Round 3)

Question

As I mentioned at the end of my Round 2 answer, I've needed to expand my code in order to produce faithfully the data needed for Figure 1 of this paper.

Unfortunately, the updates have made my script extremely slow, to the point that it doesn't finish within an amount of time I'm willing to wait.

Why is the script so slow now? This requires some explanation:

Explanation of How the Code has Changed Since Round 2

My code simulates a network of 128 leaky integrate-and-fire neurons for 10 seconds. If you look at the code posted with my Round 1 and Round 2 questions, you'll see that in those rounds I kept track of only the most recent spike (action potential) of each neuron, in a variable called t_spike_array. t_spike_array had length N = 128; it contained the time of the most recent spike for each neuron in the population.

However, if you read the paper, you'll see that, in the computation of I_syn, it calls for not just the most recent spike of each neuron, but "all the spikes emitted prior to time t by all the . . . neurons."

Hence, in the code you see below, I have replaced t_spike_array with t_spike_vector, a vector of 128 vectors, each of which starts empty. Now, when a given neuron spikes during the simulation, the time of the spike is added to this neuron's vector within t_spike_vector. In the computation of I_syn, all of the spikes in t_spike_vector are processed.

My Question

Is there anything about the way I'm performing this simulation that is inefficient? What can I do to increase the speed?

In answering these questions, here are some things to consider:

1. Note that each spike's contribution to I_syn is f(t_diff, tau_1, tau_2), where t_diff is the current time (t) minus the time of the spike (t_spike_vector[i_neur][i_spike]).

   f goes to 0 fairly quickly as t_diff increases. Hence, each spike's contribution to I_syn diminishes rapidly as time advances. Specifically, I can remove spikes that are over 130 ms old from t_spike_vector without any loss of precision in my final result.

   I assumed popping old spike times as the simulation progresses would give me a speed boost, but it doesn't seem to. In fact, in the script below, I've included code that immediately (at the next time step) pops spike times. However long I wait to pop old spikes, the script is very slow.

2. An obvious speed boost would come from representing f as a pair of differential equations and solving them using a numerical integration technique. (Note that I use a numerical integration technique -- the Euler method -- to solve for V, but solve f exactly.) My reason for not using numerical integration on f is that this paper's specific purpose is to evaluate numerical integration methods. Since the authors specify the Euler method for V, but make no mention of a numerical integration short-cut for f, I'm assuming they solved f exactly.

   This leads me to my final question: How likely is it that the authors were able to find a fast way to solve f exactly?

3. In the script below, I have n_x = 2. To produce the full dataset needed for Figure 1, I would need to set n_x to something like 100. This of course would add substantial processing time (hence why I'm using the smaller value as I refactor).

Confession

It's possible that my script is (silently) buggy. I know this site is not for debugging. Because my question regards what's making my code slow, not what's making the results incorrect -- indeed, I can't even get results due to the slowness -- I have seen it fit to post here.

(It's possible that a bug is what's making the code slow, but I believe that getting at the root of the speed issue qua speed issue is what's needed -- again justifying a post here.)

The Code (Remember, it's Very Slow)

#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>
#include <stdlib.h>
#include <sstream>
#include <stdexcept>
#include <typeinfo>
using Eigen::MatrixXd;
using Eigen::ArrayXd;


bool save_vector(const std::vector<double>& pdata,
                 const std::stringstream& file_path)
{
  std::ofstream os(file_path.str().c_str());
  if (!os.is_open())
    {
      std::cout << "Couldn't open " << file_path.str().c_str() << ".\n";
      return false;
    }
  os.precision(11);
  const int n_vals = pdata.size();
  for (int i = 0; i < n_vals; i++)
    {
      os << pdata[i];
      if (!(i + 1 == n_vals))
        {
          os << ",";
        }
    }
  os.close();
  return true;
}


std::string get_save_file()
{
  std::string dan_dir;
  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";
    }
  else if (stat("/home/despo", &statbuf) == 0 && S_ISDIR(statbuf.st_mode))
    {
      dan_dir = "/home/despo/dbliss";
    }
  std::string save_file = "/dopa_net/results/hansel/fig_1/";
  save_file = dan_dir + save_file;
  return save_file;
}


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


ArrayXd set_initial_V(const double tau, const double g_L, const double I_0,
                      const double theta, const double V_L, const int N,
                      const double c)
{
  const double T = -tau * log(1 - g_L / I_0 * (theta - V_L));
  ArrayXd V(N);
  for (int i = 0; i < N; i++)
    {
      V(i) = V_L + I_0 / g_L * (1 - exp(-c * (i - 1) / N * T / tau));
    }
  return V;
}


double square(const double base)
{
  return pow(base, 2);
}


double get_first_spike_safe(const std::vector<double>& spike_times)
{
  try
    {
      return spike_times.at(0);
    }
  catch (std::out_of_range &e)
    {
      return 999999999999;
    }
}


int main(int argc, char *argv[])
{

  // Declare and initialize constant parameters.
  const int n_x = 2;
  const double x_min = 0;   // uA / cm^2.
  const double x_max = 1;   // uA / cm^2.
  const double x_step = (x_max - x_min) / (n_x - 1);  // uA / cm^2.
  const double tau_1 = 3.0;  // ms.
  const double tau_2 = 1.0;  // ms.
  const int N = 128;
  const double dt_array[3] = {0.25, 0.1, 0.01};  // ms.
  const char* task_id = argv[argc - 1];
  const int task_id_int = task_id[0] - '0';
  const double dt = dt_array[task_id_int - 1];
  const double tau = 10;  // ms.
  const double g_L = 0.1;  // mS / cm^2.
  const double I_0 = 2.3;  // uA / cm^2.
  const double theta = -40;  // mV.
  const double V_L = -60;  // mV.
  const double c = 0.5;
  const double C = 1;  // uF / cm^2.
  const int sim_t = 10000;  // ms.

  // Declare variables set inside loops below.
  int i_t;
  int n_spikes;
  int i_saved;
  int i_spike;
  int i_neur;
  double t;
  double I_syn_bar;
  MatrixXd saved_V(N, sim_t / 2);
  double I_syn;
  double f_sum;
  ArrayXd V;
  ArrayXd A_N;
  double delta_N;
  MatrixXd V_squared;
  ArrayXd V_squared_time_mean;
  ArrayXd V_time_mean;
  ArrayXd V_time_mean_squared;
  double delta;
  double first_spike;
  double t_diff;
  std::vector<double>::iterator i_first_spike;
  std::vector <double> sigma_N_array(n_x);
  std::vector<double> spike_times;

  std::cout.precision(10);

  std::cout << "dt = " << dt << "\n";

  int i_I_syn_bar = 0;
  for (I_syn_bar = x_min; I_syn_bar < x_max; I_syn_bar += x_step)
    {

      std::cout << "I_syn_bar = " << I_syn_bar << "\n";

      // Initialize an empty vector for each neuron's spike times.
      std::vector<std::vector<double> > t_spike_vector(N,
                                                       std::vector<double>(0));

      // Process each time point in order.
      i_t = 0;
      i_saved = 0;
      for (t = 0; t < sim_t; t += dt)
        {

          if (t == 0)
            {

              // Initialize V.
              V = set_initial_V(tau, g_L, I_0, theta, V_L, N, c);

            }
          else
            {

              // Compute the sum of f over neurons and spikes.
              f_sum = 0;
              for (i_neur = 0; i_neur < N; i_neur++)
                {

                  spike_times = t_spike_vector[i_neur];
                  n_spikes = spike_times.size();
                  if (n_spikes > 0)
                    {
                      for (i_spike = 0; i_spike < n_spikes; i_spike++)
                        {
                          t_diff = t - spike_times[i_spike];
                          f_sum += f(t_diff, tau_1, tau_2);
                        }
                    }

                  // While you're here, clean out old spikes.
                  first_spike = get_first_spike_safe(spike_times);
                  while ((t - first_spike) > 0)
                    {
                      i_first_spike = t_spike_vector[i_neur].begin();
                      t_spike_vector[i_neur].erase(i_first_spike);
                      spike_times = t_spike_vector[i_neur];
                      first_spike = get_first_spike_safe(spike_times);
                    }

                }

              // Compute I_syn for this time point.
              I_syn = I_syn_bar / N * f_sum;

              // Compute V for this time point.
              V += dt * (-g_L * (V - V_L) + I_syn + I_0) / C;

            }

          // Check whether any neurons have spiked at this time point.
          for (i_neur = 0; i_neur < N; i_neur++)
            {
              if (V(i_neur) >= theta)
                {
                  t_spike_vector[i_neur].push_back(t);
                  V(i_neur) = V_L;
                }
            }

          // Save V every 1 ms from t = 5 s to t = 10 s.
          if (t >= 5000 && std::floor(t) == t)
            {
              saved_V.col(i_saved) = V;
              i_saved++;
            }

          // Increment the time counter.
          i_t++;

        }

      // Compute A_N.
      A_N = saved_V.colwise().mean();

      // Compute delta_N.
      delta_N = A_N.square().mean() - pow(A_N.mean(), 2);

      // Compute delta.
      V_squared = saved_V.unaryExpr(std::ptr_fun(square));
      V_squared_time_mean = V_squared.rowwise().mean();
      V_time_mean = saved_V.rowwise().mean();
      V_time_mean_squared = V_time_mean.square();
      delta = (V_squared_time_mean - V_time_mean_squared).mean();

      // Compute sigma_N.
      sigma_N_array[i_I_syn_bar] = delta_N / delta;

      i_I_syn_bar++;
    }

  const std::string save_file = get_save_file();
  std::stringstream complete_save_file;
  complete_save_file << save_file << dt << ".txt";
  save_vector(sigma_N_array, complete_save_file);
  complete_save_file.str("");
  complete_save_file.clear();
  return 0;

}

Answer 1

So here are a few more thoughts on your code:

Declare variables as local as possible

In general you want to declare variables in the most local scope that is possible. This makes it much easier to reason about the code and understand what is doing.

For instance spike_times is used only within a single loop iteration and should be declared there. Actually the spike_times within the for loop doesn't even have anything to do with the spike_times within the inner while loop.

The same applies for the loop variables, e.g.:

for (i_neur = 0; i_neur < N; i_neur++)

Should be

for (int i_neur = 0; i_neur < N; i_neur++)

Ther can be some exeptions for performance reasons, but in this case there are much better ways to improve performance.

Avoid unnecessary copies

spike_times = t_spike_vector[i_neur];

This copies the entire memory of the vector inside the vector, which is costly. Once it is clear (see above) that it is simply a shorthand for t_spike_vector[i_neur] used as a you can easily replace it with a const std::vector<double>& spike_times. Taking it further, you can also use a non-const reference std::vector<double>& spike_times as a shorthand even for when you write to the vector.

In C++11 this loop would be much much leaner with an

for (auto& spkike_times : t_spike_vector)

The pure C++03 way to drop the index variable are iterators:

for(std::vector<std::vector<double>>::iterator spike_it = t_spike_vector.begin(); it != t_spike_vector.end(); ++spike_it) {
    int nspikes = spike_it->size();
    ...
}

Cleaning out the old spikes.

The way you repeatedly remove the first spike is very inefficient, because std::vector::erase is a costly operation. There are two ways to approach this:

1) Use a double ended queue, std::deque. It provides efficient pop_front.

2) If you stick with the vector, at least find the right element first, and then delete all first "n" elements.

auto first_recent = std::find(spike_times.begin(), spike_times.end(), [t](auto spike){ return t <= spike; });
spike_times.erase(spike_times.begin(), first_recent);

I leave it as an exercise to convert this modern C++ code to C++03 using std::bind/std::less. Besides being much more efficient, this code also expresses your intent much better than a while loop by using std::find.

Write for clarity

To me, your code is very difficult to read in many corners, a few more expamples:

while ((t - first_spike) > 0)

Would (t > first_spike) not be much cleaner?

  for (t = 0; t < sim_t; t += dt)
    {

      if (t == 0)

Why not initialize outside of the loop? This would make the loop much easier to understand.

              if (n_spikes > 0)
                {
                  for (i_spike = 0; i_spike < n_spikes; i_spike++)

The if is absolutely redundant here. Again replace the loop with an algorithm to more clearly express your intent to accumulate/reduce a value based on all values in your vector. Use std::accumulate, just make sure to use an the initialize-value of the right type.

In C++, you do ask for permission, not forgiveness

double get_first_spike_safe(const std::vector<double>& spike_times)
{
  try
    {
      return spike_times.at(0);
    }
  catch (std::out_of_range &e)
    {
      return 999999999999;
    }
}

This is very pythonic, but not idiomatic for C++. Instead first test, then decide what to do.

Also avoid magic numbers... if you have to resort to such hacks resort to std::numeric_limits. (As disussed earlier, you don't need this hack anyway.)

On a more higher level, you main function is quite complex. While it can be improved, you should also consider to split your algorithm into more functions with clear interfaces.

At the end of the day, this may or may not improve performance. For that i refer to my comment about performance analysis tools. Based on a very short test with perf I do suspect this code is actually limited by libm / exp

By the way, you can build perf quite easy as user, but you should use the kernel version that is installed on your system.

This leads me to my final question: How likely is it that the authors were able to find a fast way to solve f exactly?

Have you asked the authors? In CS it is common to ask for practial details about papers and discuss / share research methodology with fellow researchers.