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;
}
C++11
dependency by removing the use ofauto
. 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\$ – dbliss Oct 15 '15 at 7:34auto
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 Oct 15 '15 at 11:43