Basically I'm trying to run a simulation of neurons and this is the code that gets run in the main loop of it.
class Neuron
{
private:
double v, u, a, b, c, d, current, du, dv;
public:
double x;
vector<double> pre;
int id;
vector<int> fired;
vector<int*> post;
vector<double (*)(double)> injections;
vector<double> spikes;
Neuron(Cluster *cluster, int id)
{
v = -55;
u = 0;
this->id = id;
x = uniform();
a = cluster->a();
b = cluster->b();
c = cluster->c();
d = cluster->d();
}
void integrate(double time)
{
current = 0;
for(unsigned int i = 0;i < injections.size();i++)
{
current += injections[i](time);//program crashes here!!!
}
du = a*(b*v - u);
v += .5*(0.04*pow(v, 2) + 5*v + 140 - u + current);
v += .5*(0.04*pow(v, 2) + 5*v + 140 - u + current);
u += du;
if(v >= 30)
{
v = c;
u += d;
for(unsigned int q = 0;q < post.size();q++)
{
*post[q] = 1;
}
spikes.push_back(time);
}
}
void sum()
{
for(unsigned int q = 0;q < pre.size();q++)
{
if(fired[q])
{
v += pre[q];
fired[q] = 0;
}
}
}
};
boost::barrier *neural_sync;
void main_loop(vector<Neuron*> neurons, double time_limit, double resolution)
{
for(double t = 0;t < time_limit;t += resolution)
{
printf("simulating %f\n", t);
for(unsigned int i = 0;i < neurons.size();i++)
{
neurons[i]->integrate(t);
}
neural_sync->wait();
for(unsigned int i = 0;i < neurons.size();i++)
{
neurons[i]->sum();
}
neural_sync->wait();
}
}
I'll put the full code at the end of the post if you need to see it and/or run it yourself, you'll need boost threads to build it. The program runs fine for a while and then a segmentation fault occurs at the line I point out above(line 29). The program runs that line several thousand times just fine before the fault occurs.
What I find really perplexing is that if I change
current += injections(time);
to
injections(time);
it runs fine.
Also the pointers in the injections vector point to one of two functions.
boost::mt19937 range;
boost::normal_distribution<> values_normal(0, 1);
boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > normal(range, values_normal);
double e_i(double a)
{
return 5*normal();
}
double i_i(double a)
{
return 2*normal();
}
If I alter those functions so they return smaller numbers the program ones fine.
source code:
I ran this with g++ with "g++ -o main main.cpp -lboost_thread"
main.cpp
#include <stdio.h>
#include "Brain.h"
#include <stdlib.h>
#include <time.h>
////////////////////////////////////////////////////////////////
double e_a()
{
return .02;
}
double e_b()
{
return .2;
}
double e_c()
{
return -65 + 15*pow(uniform(), 2);
}
double e_d()
{
return 8 - 6*pow(uniform(), 2);
}
double e_s()
{
return .5*uniform();
}
double e_i(double a)
{
return 5*normal();
}
///////////////////////////////////////////////////////////
double i_a()
{
return .02 + .08*uniform();
}
double i_b()
{
return .25 - .05*uniform();
}
double i_c()
{
return -65;
}
double i_d()
{
return 2;
}
double i_s()
{
return -rand()/(double)RAND_MAX;
}
double i_i(double a)
{
return 2*normal();
}
int main()
{
time_t seconds;
time(&seconds);
range.seed((unsigned int) seconds);
vector<Cluster*> clusters;
clusters.push_back(new Cluster(800, &e_a, &e_b, &e_c, &e_d));
clusters.push_back(new Cluster(200, &i_a, &i_b, &i_c, &i_d));
clusters[0]->connect(clusters[0], &e_s);
clusters[0]->connect(clusters[1], &e_s);
clusters[1]->connect(clusters[0], &i_s);
clusters[1]->connect(clusters[1], &i_s);
clusters[0]->inject(&e_i);
clusters[1]->inject(&i_i);
simulate(clusters, 1000, 1, 1);
return 0;
}
Brain.h
#ifndef NEURON_GROUP_H_INCLUDED
#define NEURON_GROUP_H_INCLUDED
#include <vector>
#include <map>
#include <math.h>
#include <fstream>
#include <boost/thread/thread.hpp>
#include <boost/thread/barrier.hpp>
#include <boost/shared_array.hpp>
#include <boost/random.hpp>
#include <boost/bind.hpp>
using namespace std;
boost::mt19937 range;
boost::uniform_real<> values(0, 1);
boost::variate_generator<boost::mt19937&, boost::uniform_real<> > uniform(range, values);
boost::normal_distribution<> values_normal(0, 1);
boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > normal(range, values_normal);
class Cluster;
class Neuron;
boost::mutex m_mutex;
class Connection
{
public:
Cluster *target;
double density, size;
double (*weights)();
Connection(Cluster *target, double (*weights)(), double density, double size)
{
this->target = target;
this->weights = weights;
this->density = density;
this->size = size;
}
};
class Injection
{
public:
double density, size;
double (*source)(double);
Injection(double (*source)(double), double density, double size)
{
this->source = source;
this->density = density;
this->size = size;
}
};
class Cluster
{
public:
vector<Neuron*> neurons;
vector<Connection*> connections;
vector<Injection*> injections;
int neuron_num, connection_num, injection_num;
double (*a)(), (*b)(), (*c)(), (*d)();
Cluster(int neuron_num, double (*a)(),double (*b)(),double (*c)(),double (*d)())
{
this->neuron_num = neuron_num;
this->a = a;
this->b = b;
this->c = c;
this->d = d;
}
void connect(Cluster *target, double (*weights)(), double density = 1, double size = 1)
{
connections.push_back(new Connection(target, weights, density, size));
}
void inject(double (*source)(double), double density = 1, double size = 1)
{
injections.push_back(new Injection(source, density, size));
}
};
class Neuron
{
private:
double v, u, a, b, c, d, current, du, dv;
public:
double x;
vector<double> pre;
int id;
vector<int> fired;
vector<int*> post;
vector<double (*)(double)> injections;
vector<double> spikes;
Neuron(Cluster *cluster, int id)
{
v = -55;
u = 0;
this->id = id;
x = uniform();
a = cluster->a();
b = cluster->b();
c = cluster->c();
d = cluster->d();
}
void integrate(double time)
{
current = 0;
for(unsigned int i = 0;i < injections.size();i++)
{
current += injections[i](time);
}
du = a*(b*v - u);
v += .5*(0.04*pow(v, 2) + 5*v + 140 - u + current);
v += .5*(0.04*pow(v, 2) + 5*v + 140 - u + current);
u += du;
if(v >= 30)
{
v = c;
u += d;
for(unsigned int q = 0;q < post.size();q++)
{
*post[q] = 1;
}
spikes.push_back(time);
}
}
void sum()
{
for(unsigned int q = 0;q < pre.size();q++)
{
if(fired[q])
{
v += pre[q];
fired[q] = 0;
}
}
}
};
boost::barrier *neural_sync;
void main_loop(vector<Neuron*> neurons, double time_limit, double resolution)
{
for(double t = 0;t < time_limit;t += resolution)
{
printf("simulating %f\n", t);
for(unsigned int i = 0;i < neurons.size();i++)
{
neurons[i]->integrate(t);
}
neural_sync->wait();
for(unsigned int i = 0;i < neurons.size();i++)
{
neurons[i]->sum();
}
neural_sync->wait();
}
}
void simulate(vector<Cluster*> brain, double time_limit, double resolution = 1, int thread_num = 1)
{
vector<vector<Neuron*> > thread_neurons(thread_num);
Cluster *tmp;
unsigned int thread = 0;
int id = 0, total_neurons = 0;
bool cull_threads = true;
Neuron *neuron;
for(unsigned int i = 0;i < brain.size();i++)
{
tmp = brain[i];
total_neurons += tmp->neuron_num;
for(int u = 0;u < tmp->neuron_num;u++)
{
neuron = new Neuron(tmp, id++);
tmp->neurons.push_back(neuron);
thread_neurons[thread].push_back(neuron);
thread = (thread + 1) % thread_num;
if(thread == 0)
cull_threads = false;
}
}
if(cull_threads)
{
thread_num = thread;
thread_neurons.resize(thread_num);
}
printf("building neurons\n");
double target_x;
Neuron *from_neuron, *to_neuron;
for(unsigned int i = 0;i < brain.size();i++)
{
tmp = brain[i];
for(unsigned int u = 0;u < tmp->neurons.size();u++)
{
from_neuron = tmp->neurons[u];
for(unsigned int q = 0;q < tmp->connections.size();q++)
{
target_x = uniform();
for(unsigned int p = 0;p < tmp->connections[q]->target->neurons.size();p++)
{
to_neuron = tmp->connections[q]->target->neurons[p];
if(from_neuron != to_neuron && abs(target_x - to_neuron->x) <= tmp->connections[q]->size && uniform() < tmp->connections[q]->density)
{
to_neuron->pre.push_back(tmp->connections[q]->weights());
to_neuron->fired.push_back(0);
from_neuron->post.push_back(&to_neuron->fired[to_neuron->fired.size() - 1]);
}
}
}
for(unsigned int q = 0;q < tmp->injections.size();q++)
{
target_x = uniform();
if(abs(target_x - from_neuron->x) <= tmp->injections[q]->size && uniform() < tmp->injections[q]->density)
{
from_neuron->injections.push_back(tmp->injections[q]->source);
}
}
}
}
neural_sync = new boost::barrier(thread_num);
printf("starting simulation\n");
boost::shared_array<boost::thread> threads = boost::shared_array<boost::thread>(new boost::thread[thread_num]);
for(int i = 0;i < thread_num;i++)
{
threads[i] = boost::thread(boost::bind(&main_loop, thread_neurons[i], time_limit, resolution));
}
for(int i = 0;i < thread_num;i++)
{
threads[i].join();
}
printf("finished simulation\n");
delete neural_sync;
printf("building results\n");
/*fstream file_x("trace_x.bin", ios::out|ios::binary);
fstream file_y("trace_y.bin", ios::out|ios::binary);
for(unsigned int i = 0;i < brain.size();i++)
{
tmp = brain[i];
for(unsigned int u = 0;u < brain[i]->neurons.size();u++)
{
for(unsigned int q = 0;q < brain[i]->neurons[u]->spikes.size();q++)
{
file_y.write((char*)&brain[i]->neurons[u]->id, sizeof(int));
file_x.write((char*)&brain[i]->neurons[u]->spikes[q], sizeof(double));
}
}
}
file_x.close();
file_y.close();*/
printf("built results\n");
}
#endif