dneuron.cpp

a useful spiking neural networks simulator

}

template <bool debug>
bool tneuron_ext<debug>::broadcastMessage(ntime_t t, message_base *mess)
{
	bool ret = false;
	for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		ret = (*i)->processMessage(t, *mess) || ret;
	return ret;
}

template <bool debug>
real tneuron_ext<debug>::calcSignal( ntime_t current_time )
{
	real signal = 0.0;
	for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
	{
		signal += (*i)->getValue(current_time);
	}
	return signal;
}

template <bool debug>
void tneuron_ext<debug>::fire(tscheduler &scheduler, ntime_t current_time)
{
	{
		prof<1> profobj("Fire");
		totalfire++;

		for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
		{
	//#ifdef DEBUG
	//		cout << "schedule pulse at " << current_time + i->getDelay() << " to " << i->getDest().getName() << endl;
	//#endif
			scheduler.scheduleEvent(current_time + (*i)->getDelay(), **i);
		}
		last_fire = current_time;
		for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		{
			(*i)->setZeroPoint(current_time);
		}
	}
	scheduleFire(scheduler, current_time);
}


struct valdomain { real value, upslope, ceil, downslope, floor; };

ostream& operator<<(ostream &o, const valdomain &v) 
{
	return o << setprecision(8) << "val=" << v.value << ", u=" << v.upslope << ", c=" << v.ceil << ", d=" << v.downslope << ", f=" << v.floor;
}

const bool debug_domain = false;

real domain_floor_time( vector<valdomain> &sig_domain, ntime_t next_incontinuity )
{
	if( debug_domain )
		cout << "enter domain_floor" << endl;
	ntime_t search_from = 0;
	ntime_t old_search_from;

	do {
		if( debug_domain )
			cout << "domain_floor search_from=" << search_from << endl;

		real sig = 0.0, slo = 0.0;
		for( vector<valdomain>::iterator i = sig_domain.begin(); i != sig_domain.end(); i++ )
		{
			real nsig = i->value + i->downslope * search_from;
			if( debug_domain )
				cout << "domain nsig=" << nsig << ", " << *i << endl;
			if( nsig < i->floor ) {
				sig += i->floor;
			}
			else {
				sig += nsig; 
				slo += i->downslope;
			}
		}
		old_search_from = search_from;
		if( debug_domain )
			cout << "domain_floor sig=" << sig << ", slo=" << slo << endl;
		if( fabs(slo) < 1e-8 )
		{
			if( fabs(sig) > 1e-4 )
				search_from = Infinity;
		}
		else
			search_from -= sig / slo;
	} while( old_search_from < search_from && search_from < next_incontinuity);

	if( debug_domain )
		cout << "leave domain_floor search_from=" << search_from << endl;
	return search_from;
}

real domain_ceil_time( vector<valdomain> &sig_domain, ntime_t next_incontinuity )
{
	if( debug_domain )
		cout << "enter domain_ceil" << endl;
	ntime_t search_from = 0;
	ntime_t old_search_from;

	do {
		if( debug_domain )
			cout << "domain_ceil search_from=" << search_from << endl;
		real sig = 0.0, slo = 0.0;
		for( vector<valdomain>::iterator i = sig_domain.begin(); i != sig_domain.end(); i++ )
		{
			if( debug_domain )
				cout << "domain " << *i << endl;
			real nsig = i->value + i->upslope * search_from;
			if( nsig >= i->ceil ) {
				sig += i->ceil;
			}
			else {
				sig += nsig; 
				slo += i->upslope;
			}
		}
		old_search_from = search_from;
		if( debug_domain )
			cout << "domain_ceil sig=" << sig << ", slo=" << slo << endl;
		if( fabs(slo) < 1e-8 )
		{
			if( fabs(sig) > 1e-4 )
				search_from = Infinity;
		}
		else
			search_from -= sig / slo;
	} while( old_search_from < search_from && search_from < next_incontinuity);

	if( debug_domain )
		cout << "leave domain_ceil search_from=" << search_from << endl;
	return search_from;
}

template <bool debug>
void tneuron_ext<debug>::pulseArrive(tscheduler &scheduler, ntime_t current_time, real pulse_level)
{
	real curval = pulses->getValue(current_time);

	pulses->setParam(curval + pulse_level, current_time);

	if( getDeb() )
	cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
		<< " Pulse arrived, signal level " << pulse_level << ", new pulse " << pulses->getValue(current_time) << endl; 

	scheduleFire(scheduler, current_time);
}

template <bool debug>
void tneuron_ext<debug>::pulseArrive(tscheduler &scheduler, ntime_t current_time, message_base *base)
{
//	real curval = pulses->getValue(current_time);

	if( getDeb() )
	cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
		<< " Pulse arrived, message=" << base->getMessageId() << flush;
	bool ret = sendMessage(current_time, base);
	if( getDeb() )
		cout << " result=" << ret << endl;

	scheduleFire(scheduler, current_time);
}

template <bool debug>
void tneuron_ext<debug>::pulseArrive(tscheduler &scheduler, ntime_t current_time, func_base *base)
{
//	real curval = pulses->getValue(current_time);

	if( getDeb() )
	cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
		<< " Pulse arrived, func=" << base->getDescription() << flush;
	addExt(base);

	scheduleFire(scheduler, current_time);
}

template <bool debug>
void tneuron_ext<debug>::scheduleFire(tscheduler &scheduler, ntime_t current_time, bool /*resched*/)
{
	prof<1> profobj("scheduleFire");
//	if( getDeb() ) cout << "enter tneuron_ext::schedulefire" << endl;
	
	real signal = 0.0;
	real deriv1st = 0.0;
	real deriv2nd = 0.0;
	real next_incontinuity = Infinity;
	ntime_t next_known_pulse;

	int extssize = exts.size();
	vector<valdomain> sig_domain(extssize);
	vector<valdomain> d1st_domain(extssize);
	vector<valdomain> d2nd_domain(extssize);

	last_simulate = current_time;
	real maxgrad = 0.0;

	if( use_next_known_pulse ) { ////////////////////////////////////////////
		prof<2> profobj("tneuron_ext::0 calcval");
		while( !queue()->empty() && &(queue()->top().getAction()) == this )
		{
			totalrescheduled++;
			queue()->pop();
		}
		if( queue()->empty() )
		 	next_known_pulse = Infinity;
		else
			next_known_pulse = queue()->top().getTime();

		vector<valdomain>::iterator sdi = sig_domain.begin();
		for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		{
			real v = (*i)->getValue(current_time);
			sdi++->value = v;
			signal += v;

//			cout << "v=" << v << ", m=" << (*i)->getMaxGradient(current_time) << ", desc= "<< (*i)->getDescription() << endl;

			maxgrad += (*i)->getMaxGradient(current_time);
		}
	} ///////////////////////////////////////////////////////////////////////
	else {
		next_known_pulse = Infinity;
		vector<valdomain>::iterator sdi = sig_domain.begin();
		for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		{
			real v = (*i)->getValue(current_time);
			valdomain vd; vd.value = v;
			sdi++->value = v;
			signal += v;
	//cout << "v=" << v << ", desc= "<< (*i)->getDescription() << endl;
		}
	}

	if( signal >= - 1e-8 )
	{
		if( getDeb() )
			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
				<< " Fire, signal level " << signal << endl; 
		fire(scheduler, current_time);
	}
	else
	{
		if( use_next_known_pulse ) { ////////////////////////////////////////////
			if( signal < - maxgrad * (next_known_pulse - current_time) )
			{
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " Stay, signal level " << signal << ", filtered by max gradient " << maxgrad << ", next " << next_known_pulse << endl;
				totalfiltered_maxgrad++;
				return;
			}
		}
		{
		prof<2> profobj("tneuron_ext::1 sigdomain");
		vector<valdomain>::iterator sdi = sig_domain.begin();
		for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		{
	//		if( getDeb() )
	//			cout << "func_base " << (*i)->getDescription() << endl;
			while( (*i)->shouldDelete(current_time) )
			{
				prof<2> profobj("delete");
				delete *i;
				exts.erase(i);
				sig_domain.erase(sdi);
				if( i == exts.end() )
					goto EXIT_DOM_LOOP;
			}

	//		if( sdi == sig_domain.end() ) cout << "Why? " << sig_domain.size() << " / " << exts.size() << endl;
			(*i)->getValueDomain(current_time, sdi->upslope, sdi->ceil, sdi->downslope, sdi->floor);
			sdi++;
		}
		EXIT_DOM_LOOP: ;
		}
		real value_bound;
		{
		prof<2> profobj("tneuron_ext::2 filter");
		value_bound = domain_ceil_time(sig_domain, next_incontinuity - current_time);

		if( value_bound + current_time >= next_incontinuity && next_incontinuity < Infinity )
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Stay, signal level " << signal << ", Reschedule at next incontinuity " << next_incontinuity << endl; 
			setLoopBack(scheduler, next_incontinuity);
			totalfiltered_incontinuity++;
			return;
		}
		if( use_next_known_pulse && value_bound + current_time >= next_known_pulse && next_known_pulse < Infinity )
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Stay, signal level " << signal << ", wait for the next known pulse at " << next_known_pulse << endl; 
			totalfiltered_nextpulse++;
			return;
		}
		}
		real next_bound, bound_val;
		{
		prof<2> profobj("tneuron_ext::3 bounder");
		if( iporder >= 1 )
		{
			vector<valdomain>::iterator d1d = d1st_domain.begin();
			for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
			{
				real v = (*i)->get1stDeriv(current_time);