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

I realize this is quite an odd/complicated problem, but can you reduce the amount of code necessary to produce it? Maybe you can play around and remove as many lines as possible and still produce the error. Often that will help you isolate the bug yourself, or if not, then it's much easier for us to dig in.

I don't know if this is drastically simpler, I took almost 70 lines out of brain.h but it is now strictly a single threaded program.

I wanted to make it simpler but I found two simplifications that made the program work, see below.

Brain.h

#ifndef NEURON_GROUP_H_INCLUDED
#define NEURON_GROUP_H_INCLUDED

#include <vector>
#include <map>
#include <math.h>

#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;

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
{
	public:
		double v, u, a, b, c, d, current, du, dv;
		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;
				}
			}
		}
};

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);
		}
		for(unsigned int i = 0;i < neurons.size();i++)
		{
			neurons[i]->sum();
		}
	}
}

void simulate(vector<Cluster*> brain, double time_limit, double resolution = 1)
{
    vector<Neuron*> neurons;
	Cluster *tmp;
	int id = 0, total_neurons = 0;
	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);
			neurons.push_back(neuron);
		}
	}

printf("building neurons\n");
	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++)
			{
				for(unsigned int p = 0;p < tmp->connections[q]->target->neurons.size();p++)
				{
					to_neuron = tmp->connections[q]->target->neurons[p];
                    to_neuron->pre.push_back(tmp->connections[q]->weights());
                    to_neuron->fired.push_back(0);
                    from_neuron->post.push_back(&to_neuron->fired.back());
				}
			}
			for(unsigned int q = 0;q < tmp->injections.size();q++)
			{
                from_neuron->injections.push_back(tmp->injections[q]->source);
			}
		}
	}

printf("starting simulation\n");

	main_loop(neurons, time_limit, resolution);

printf("finished simulation\n");
}

#endif

if I removed lines 107-120 or if I remove lines 178-187 the program runs fine.

I've gotten it to work now. I'm really not sure what was wrong with what I did but I also have to admit I was treading programming waters I didn't under stand too well.

I wanted to have a multithreaded program so every neuron in my simulation had a vector of integers, one integer for each neuron that could send input to it, this vector was called fired. Each neuron also had a vector of integer pointers, and each pointer in these vectors pointed to an integer in the fired vector of each neuron they could send input to, call this vector post. So when a neuron sent input to another neuron it dereferenced each pointer in in the post vector and set the integer to 1. Next, each neuron would check the integers in the fired vector to see if a neuron fired and take appropriate action.

Strangely, I don't see how any of this could be related to my original problem...but the program works.

This question has already been answered. Start a new discussion instead.