dneuron.cpp

a useful spiking neural networks simulator

/** 
 * @file  dneuron.cpp
 * @brief Neurons
 *
 * @author Makino, Takaki <[e-mail]>
 * @date 2003-05-01
 * @version $Id: dneuron.cpp,v 1.7 2003/05/07 09:32:23 t Exp $
 *
 *  Copyright (C) 2003 Makino, Takaki. 
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2, or (at your option)
 *  any later version.

 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  
 */

#include "dneuron.h"
#include <math.h>
#include <mak/profile.h>

using namespace mak;

//#define DEBUG

/////////////////////////////////////////////////////////////////////////
namespace punnets_common {
/////////////////////////////////////////////////////////////////////////

u_int64_t totalfire = 0;
u_int64_t totalpulse = 0;
u_int64_t totalpartition = 0;
u_int64_t totalpartition_nonewton[4] = {0,0,0,0};
u_int64_t totalpartition_newton[4] = {0,0,0,0};
u_int64_t totalpeaksearch[3] = {0,0,0};
u_int64_t totalpeakenclosing = 0;
u_int64_t totalrescheduled = {0};
u_int64_t totalfiltered_maxgrad = 0;
u_int64_t totalfiltered_incontinuity = 0;
u_int64_t totalfiltered_nextpulse = 0;

const int iporder = 2;
const ntime_t ip_delta_t = 1e-5;

//const bool use_next_known_pulse = false;
const bool use_next_known_pulse = true;

using namespace std;

class tpulse : public taction
{
	tneuron_base &dest;
	real level;

public:
	tpulse(tneuron_base &idest, real ilevel) : taction(), dest(idest), level(ilevel) { }
	
	virtual void activate(tscheduler &scheduler, ntime_t current_time) 
		{ 
//		if( getDeb() )
//			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << dest.getName() << " Pulse Arrived (external input)" << endl;
		dest.pulseArrive( scheduler, current_time, level ); delete this; }

	tqueue *queue() const { return dest.queue(); }
	virtual const char *getClassName() const { return "tpulse"; };
};

class tmessage : public taction
{
	tneuron_base &dest;
	message_base *mess;
	real level;

public:
	tmessage(tneuron_base &idest, message_base *imess) : taction(), dest(idest), mess(imess) { }
	virtual ~tmessage() { delete mess; }
	
	virtual void activate(tscheduler &scheduler, ntime_t current_time) 
		{ 
//		if( getDeb() )
//			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << dest.getName() << " Pulse Arrived (external message " << mess->getMessageId() << ")" << endl;
		dest.pulseArrive( scheduler, current_time, mess ); delete this; }

	tqueue *queue() const { return dest.queue(); }
	virtual const char *getClassName() const { return "tmessage"; };
};


taction &makePulse( tneuron_base &idest, real ilevel )
{
	return *new tpulse( idest, ilevel );
}

taction &makePulse( tneuron_base &idest, message_base *mess )
{
	return *new tmessage( idest, mess );
}


}  // namespace punnets

namespace punnets_private {

template <bool debug>
tneuron<debug>::~tneuron()

{
	for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
		delete *i;
}


template <bool debug>
tneuron<debug>::tneuron(string iname /*= ""*/, 
				 ntime_t isig_hv_period /*= def_sig_hv_period */,
				 ntime_t ithr_hv_period /*= def_thr_hv_period */,
				 real imin_threshold    /*= def_min_threshold */,
				 real imax_threshold    /*= def_max_threshold */ )
		: tneuron_base(iname), 
		  sig_level(0.0), last_simulate(0.0), last_fire(0.0),
		  sig_hv_period(isig_hv_period), thr_hv_period(ithr_hv_period),
		  min_threshold(imin_threshold), max_threshold(imax_threshold),
		  coeff_sigdecay(- M_LN2 / isig_hv_period), 
		  coeff_thrdecay( - M_LN2 / ithr_hv_period )
{ }


template <bool debug>
void tneuron<debug>::fire(tscheduler &scheduler, ntime_t current_time)
{
	prof<1> pf("tneuron::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);
	}
	sig_level = 0.0;
	last_fire = current_time;
}

template <bool debug>
void tneuron<debug>::scheduleFire(tscheduler &scheduler, ntime_t current_time, bool /*resched*/)
{
	prof<1> pf("tneuron::scheduleFire");
	real current_threshold = 0.0;
	if( sig_level >= min_threshold )
	{
		// 嵟戝偟偒偄抣傪挻偊偰偄偨傜懄嵗偵敪壩
		if( sig_level >= max_threshold || (current_threshold = getCurrentThrLevel( current_time )) - sig_level < epsilon )
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Fire, level " << sig_level << ", threshold " << current_threshold << endl;
			fire(scheduler, current_time);
			return;
		}
		else
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Candid, level " << sig_level << ", threshold " << current_threshold << endl;
		}
	}
	else
		if( getDeb() )
			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Accumulate, " << sig_level << endl;
}

/************************************************************************/

template <bool debug>
void tneuron_ext_const<debug>::scheduleFire(tscheduler &scheduler, ntime_t current_time, bool resched)
{
	real current_threshold = 0.0;
	if( sig_level >= min_threshold || sig_converge_level > min_threshold )
	{
		// 嵟戝偟偒偄抣傪挻偊偰偄偨傜懄嵗偵敪壩
		if( sig_level >= max_threshold || (current_threshold = getCurrentThrLevel( current_time )) - sig_level < epsilon )
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Fire, level " << sig_level << ", threshold " << current_threshold << endl;
			fire(scheduler, current_time);
			if( sig_converge_level < min_threshold )
				return;
			current_threshold = max_threshold;
		}
		else
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Candid, level " << sig_level << ", threshold " << current_threshold << endl;
			if( ! resched )
				return;
		}

		// 崱屻丄偟偒偄抣埲壓偱偁傟偽丄挻偊傞壜擻惈偑偁傞偐挷傋傞
		real sig_decay_level = sig_level - sig_converge_level;
		real thr_decay_level = current_threshold - min_threshold;
		// 嵟彫偟偒偄抣偑怣崋廂懇抣枹枮側傜妋幚偵墇偊傞偺偱丄偦傟埲奜偺応崌偩偗挷傋傞
		if( min_threshold >= sig_converge_level )
		{
			// 怣崋偺敿尭婜偺傎偆偑挿偄応崌埲奜偼壜擻惈偑側偄
			if( coeff_sigdecay - coeff_thrdecay <= 0 )
			{
				if( getDeb() )
					cout << current_time << " Fail neuron " << name << endl;
				return;
			}
				
			// (偟偒偄抣娭悢 - 怣崋娭悢) 偑嬌戝偵側傞帪崗
			ntime_t maximal_time = 
				(coeff_sigdecay - coeff_thrdecay) *
				log( ( coeff_thrdecay * thr_decay_level ) / 
				     ( coeff_sigdecay * sig_decay_level )   );
			// 嬌戝偵側偭偨偲偒偵墇偊偰偄傞偐丠
			real siglvl_at_maximal = sig_decay_level * exp(coeff_sigdecay * maximal_time) + sig_converge_level;
			real thrlvl_at_maximal = thr_decay_level * exp(coeff_thrdecay * maximal_time) + min_threshold;
			if( siglvl_at_maximal < thrlvl_at_maximal )
			{
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Fail, max_lvl " << siglvl_at_maximal << ", threshold " << thrlvl_at_maximal << endl;
				return;	// 墇偊偰偄側偄丄敪壩偺壜擻惈側偟
			}
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Possible, max_lvl " << siglvl_at_maximal << ", threshold " << thrlvl_at_maximal << endl;
//			cout << T << endl;
		}
		else
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Possible" << endl;
//			cout << T << endl;
		}
		// 偟偽傜偔屻偱敪壩偡傞偺偱偦偺帪崗 t 傪挷傋傞
		// 僯儏乕僩儞朄 T = exp(t)

		// 弶婜抣偼 1.0 偩偑丄侾夞傑傢偟偨屻偺抣偱寁嶼
		real T = 1.0 - (sig_level - current_threshold) / (coeff_sigdecay * sig_decay_level - coeff_thrdecay * thr_decay_level);

		int count = 0;

		ntime_t fire_time = log(T);
		while(1) {
			real atb = sig_decay_level * exp(fire_time * coeff_sigdecay);
			real ctd = thr_decay_level * exp(fire_time * coeff_thrdecay);
			real T_ = (atb - ctd - min_threshold + sig_converge_level) / (coeff_sigdecay * atb - coeff_thrdecay * ctd);
			if( fabs(T_ ) < epsilon )
				break;
			T *= 1.0 - T_;
//				cout << T << "\t" << T_ << endl;
			fire_time = log(T);
			count++;
		}

		// 敪壩帪崗偵傕偆堦搙僗働僕儏乕儕儞僌
		if( getDeb() )
		{
			real siglvl_at_fire = sig_decay_level * exp(coeff_sigdecay * fire_time) + sig_converge_level;
			real thrlvl_at_fire = thr_decay_level * exp(coeff_thrdecay * fire_time) + min_threshold;
			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Reschedule, time " << fire_time << ", fire_lvl " << siglvl_at_fire << ", threshold " << thrlvl_at_fire << ", newton " << count << endl;
		}
//				cout << "newton " << count << endl;
		scheduler.scheduleEvent(current_time + fire_time, *this);
	}
	else
		if( getDeb() )
			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Accumulate, level " << sig_level << endl;
}

/************************************************************************/


/************************************************************************
 ***** tneuron_ext ******************************************************
 ************************************************************************/

template <bool debug>
tneuron_ext<debug>::~tneuron_ext()
{
	for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
		delete *i;
	for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		delete *i;
}

template <bool debug>
tneuron_ext<debug>::tneuron_ext(string iname /*= ""*/, real sig_hv_period)
		: tneuron_base(iname), 
		  last_simulate(-10000000), last_fire(-10000000),
		  lambda(M_LN2 / sig_hv_period) 
{
	pulses = new func_delta_int(0, 0);
	addExt(pulses);
}

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