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

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?

• switching to a vector got me to a total computing time of ~5s as compared to the original ~7s. (this is still far too slow.) after this change, i'm spending ~88.8% of the time in get_f. i'm wondering whether i need a larger algorithmic change to get_I_syn. maybe i should pre-compute values for get_f for the full range of values for t-spike_times[i] i'm expecting to get? – dbliss Nov 21 '15 at 3:54
• for reference, the paper i'm replicating, from 1998 (see the Round 1 question for a link to the paper), reports that it takes ~5s to run a 10-s simulation of 138 neurons. it takes me this long to simulate just 1 neuron for 10 s. – dbliss Nov 21 '15 at 3:59
• i should also mention that 75.3% of the time is spent in exp -- so it's really the use of this math function over and over again that's slowing me down. – dbliss Nov 21 '15 at 4:14
• is this the full working example that takes 7s? What compiler with what options are you using on what kind of a machine? On my i7-4790K with gcc5.2 -O3 it runsin 0.25s. – Zulan Nov 21 '15 at 9:34
• @Zulan yep, this is the full working example. i get 8 s using the -O3 compiler option. i'm not sure exactly what to tell you about my machine. this is the processor info: Intel® Core™2 Duo CPU T5250 @ 1.50GHz × 2. it's an old machine, but not 1998 old. – dbliss Nov 22 '15 at 20:58

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;
}
}
}
• thanks very much. your modification of get_f is exactly the sort of thing i wanted to do but wasn't sure how to achieve. regarding get_I_syn, yes, i was aware that it always returns 0. the script i posted, though complete in a sense, is really a toy version of the full script i'm going for, which will compute get_I_syn for a wide range of values of I_syn_bar. (your assumptions about get_f, though, hold for the full script, so this answer is very, very helpful.) – dbliss Nov 22 '15 at 21:54
• the code in my question uses the Euler method for numerical integration. i've supplemented this code with an "exact" method, where, rather than go forward in time in steps of Delta_t, i analytically solve for the next spike time in the network, and move straight to that time. this means that the spike times are no longer multiples of Delta_t. hence, it would seem a cache of f values no longer helps me. what would your advice be for this situation? using exp is still very slow. i'm currently testing "fast" replacements for exp. but maybe the fastest thing would be to vectorize? – dbliss Dec 5 '15 at 0:04
• there are also probably compiler options that would speed things up that i haven't considered. my google searches have revealed that my problem isn't unique, but i'm not seeing a canonical solution anywhere. it'd be nice to have one. i'd rather not, but maybe i should post a Round 5 . . . – dbliss Dec 5 '15 at 0:13