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

Question

You may want to take a look at Rounds 1, 2, and 3, though that isn't necessary for understanding what's below.

The major change since Round 3 is that my code is much cleaner and I'm including profiling information here.

Here's my script:

#include "hansel.h"
#include <deque>


int main()
{

  // Parameters.
  const double tau = 10;
  const double g_l = 0.1;
  const double I_0 = 2.3;
  const double theta = -40;
  const double V_l = -60;
  const double c = 0.5;
  const double N = 128;
  const double I_syn_bar = 0;
  const double tau_1 = 3;
  const double tau_2 = 1;
  const double C = 1;
  const double Delta_t = 0.25;

  // Simulation.
  std::deque<double> spike_times;
  double V = set_initial_V(tau, g_l, I_0, theta, V_l, c, N);
  for (double t = 0; t < 10000; t += Delta_t)
    {
      double I_syn = get_I_syn(I_syn_bar, N, tau_1, tau_2, t, spike_times);
      double dV_dt = get_dV_dt(g_l, V, V_l, I_syn, I_0, C);
      V += Delta_t * dV_dt;
      if (V > theta)
        {
          V = V_l;
          spike_times.push_back(t + Delta_t);
        }
    }

}

Here's the header file it includes:

#include <math.h>
#include <deque>


const double get_T(const double tau, const double g_l, const double I_0,
                   const double theta, const double V_l)
{
  return -tau * log(1 - g_l / I_0 * (theta - V_l));
}


double set_initial_V(const double tau, const double g_l, const double I_0,
                     const double theta, const double V_l, const double c,
                     const double N)
{
  const double T = get_T(tau, g_l, I_0, theta, V_l);
  const double i = 1;
  return V_l + I_0 / g_l * (1 - exp(-c * (i - 1) / N * T / tau));
}


double get_dV_dt(const double g_l, const double V, const double V_l,
                 const double I_syn, const double I_0, const double C)
{
  return (-g_l * (V - V_l) + I_syn + I_0) / C;
}


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


double get_I_syn(const double I_syn_bar, const double N, const double tau_1,
                 const double tau_2, const double t,
                 const std::deque<double> & spike_times)
{
  double sum_f = 0;
  const int n_spikes = spike_times.size();
  for (int i = 0; i < n_spikes; i += 1)
    {
      sum_f += get_f(tau_1, tau_2, t - spike_times[i]);
    }
  return I_syn_bar / N * sum_f;
}

Profiling the script, I see that ~100% of the time is spent in calls to get_I_syn. Within that, time is split between calls to get_f (50.1%) and use of the deque [] operator (45.4%).

I'm hoping to speed up this script by a factor of ~10.

My questions are

1) Is there a way to refactor this line to increase speed:

return 1 / (tau_1 - tau_2) * (exp(-t / tau_1) - exp(-t / tau_2));

2) Is there a way to vectorize (or otherwise speed up) this loop:

for (int i = 0; i < n_spikes; i += 1)
  {
     sum_f += get_f(tau_1, tau_2, t - spike_times[i]);
  }

3) Is it likely to be faster to use an array or a vector instead of the deque object?

Answer 1

get_I_syn() currently always returns 0

I'm not sure if you are aware of this, but because I_syn_bar is a constant 0, this makes get_I_syn() always return 0. This means you can just get rid of all calls to get_I_syn(), which is where 100% of your time is being spent.

However, assuming that there is some mixup here and I_syn_bar is supposed to be non-zero, I will explain how you can speed your program up by a factor of 50+.

Caching return values of get_f()

If you look at the get_f() function, everything is a constant except for t, which is passed in as t - spike_times[i] from get_I_syn(). If you look closely at what t and spike_times[i] are, they are all multiples of Delta_t. In fact, t - spike_times[i] can only take on one of 10000 values, which are the first 10000 multiples of Delta_t.

Knowing this, you can precompute the 10000 possible values of get_f(). Then inside get_I_syn(), you can look up the precomputed values instead of doing the expensive computation.

I modified your program to do this and the program sped up by more than 50x.

Modified program

Here is the modified program. I copied the header contents directly into the c++ file for my own convenience. As you can see, I increased the number of iterations from 10000 to 300000 because otherwise the program ran too fast.

#include <deque>
#include <math.h>

const double get_T(const double tau, const double g_l, const double I_0,
                   const double theta, const double V_l)
{
  return -tau * log(1 - g_l / I_0 * (theta - V_l));
}


double set_initial_V(const double tau, const double g_l, const double I_0,
                     const double theta, const double V_l, const double c,
                     const double N)
{
  const double T = get_T(tau, g_l, I_0, theta, V_l);
  const double i = 1;
  return V_l + I_0 / g_l * (1 - exp(-c * (i - 1) / N * T / tau));
}


double get_dV_dt(const double g_l, const double V, const double V_l,
                 const double I_syn, const double I_0, const double C)
{
  return (-g_l * (V - V_l) + I_syn + I_0) / C;
}


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


double get_I_syn(const double I_syn_bar, const double N, const double *f_cache,
                 const int t_index, const int *spike_times, int n_spikes)
{
    double sum_f = 0;
    for (int i = 0; i < n_spikes; i++)
        sum_f += f_cache[t_index - spike_times[i]];
    return I_syn_bar / N * sum_f;
}

#define ITERATIONS        300000

// Declare in global scope to not overflow the stack.
static double f_cache[ITERATIONS];
static int spike_times[ITERATIONS];

int main()
{
    // Parameters.
    const double tau = 10;
    const double g_l = 0.1;
    const double I_0 = 2.3;
    const double theta = -40;
    const double V_l = -60;
    const double c = 0.5;
    const double N = 128;
    const double I_syn_bar = 2;
    const double tau_1 = 3;
    const double tau_2 = 1;
    const double C = 1;
    const double Delta_t = 0.25;
    int    i = 0;
    double t = 0;

    // Simulation.
    int numSpikes = 0;
    double V = set_initial_V(tau, g_l, I_0, theta, V_l, c, N);

    // Precompute all the possible values returned by get_f().
    for (i = 0, t = 0; i < ITERATIONS; i++, t += Delta_t)
        f_cache[i] = get_f(tau_1, tau_2, t);

    for (i = 0; i < ITERATIONS; i++)
    {
        double I_syn = get_I_syn(I_syn_bar, N, f_cache, i, spike_times,
                numSpikes);
        double dV_dt = get_dV_dt(g_l, V, V_l, I_syn, I_0, C);
        V += Delta_t * dV_dt;
        if (V > theta)
        {
            V = V_l;
            spike_times[numSpikes++] = i + 1;
        }
    }
}