dtest.cpp

a useful spiking neural networks simulator

/** 
 * @file  dtest.cpp
 * @brief Test program of the punnets library.
 *
 *   This program tests the performance of delayed firing simulations.
 *   Every neuron has an sinusoidal external input and interconnecting synapses, 
 *   which causes alpha-function-style (difference of two exponential functions)
 *   response.  Two hundred pulses are injected to produce initial activity.
 *   You can change parameters by command-line options; try <code>dtest -h</code>
 *   to show the list of options.
 *
 * @author Makino, Takaki <[e-mail]>
 * @date 2003-05-01
 * @version $Id: dtest.cpp,v 1.7 2003/05/08 08:23:56 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 <iostream>
#include <fstream>
#include <sstream>
#include "dneuron.h"
#include "dsched.h"
#include "dlogger.h"

#include <mak/cmdopt.h>

using namespace std;
using namespace punnets_nodebug;
using namespace mak;

/// An observer neuron class.
///
/// When this neuron received a pulse, it counts the number of pulses.
/// If constructed with true <var>show</var> parameter, it also reports the arrival of a pulse to the console.

class tobserver : public tneuron_base
{
	int npulses;
	bool show;
public:
	/// The constructor.
	tobserver(string iname, bool ishow = false) : tneuron_base(iname), show(ishow) { npulses = 0; }
	/// When a pulse is arrived, it counts and reports the pulse arrival.
	virtual void pulseArrive(tscheduler &, ntime_t current_time, real pulse_level)
	{
		npulses++;
		if( show )
			cout << current_time << " '" << name << "' Pulse observed " << pulse_level << "  fire=" << totalfire << ", pulse=" << totalpulse << endl;
	}

	/// Get the number of pulses observed.
	int getNPulses() const { return npulses;}
	/// Clear the pulse counter.
	void clearNPulses() { npulses = 0; }
};


int main(int argc, char **argv)
{
	cout << "Punnets performance test program  (c) 2003 Makino Takaki." << endl;
	cout << "Punnets is free software, covered by the GNU General Public License." << endl;
	cout << endl;

	cmdopt_parser p;
	cmdopt<int> nneurons("The number of neurons in the simulation.",
			'n', "neurons", 100, 0, 1000);
	cmdopt<int> nlogging("The number of neurons to be logged. (Default 5)",
			'l', "log", 5, 0, 1000);
	cmdopt<double> stepsize("Specify the time step for logging. (Default 0.125)",
			's', "step", 0.125, 0, 1000);
	p.add( New<cmdopt_help>("Show this help.", 
			'h', "help") );
	p.add(Ref(nneurons)); p.add(Ref(nlogging)); 
	p.parse(argc, argv);
	
	const unsigned int nconnections = 10;
	const ntime_t hv_time = 5;
	const ntime_t simulate_until = 1000;
	const char * const datafilename = "test.dat";
	const char * const plotfilename = "test.plt";

	cout << "Preparing network..." << endl;

	vector<tneuron_ext *> exts(nneurons.get());

	tobserver obs("ob");

	srand48(12345);

	tscheduler scheduler;

	cout << "Produce simulation log file '" << datafilename << "'." << endl;
	ofstream ofs(datafilename);
	tlogger logger(ofs, stepsize.get(), stepsize.get(), simulate_until);
	logger.schedule(scheduler);
	message_base *mess = new func_exp_diff::message_add_event_time;

	for( int i=0; i<nneurons.get(); i++ )
	{
		ostringstream oss;
		oss << "n" << i;
		exts[i] = new tneuron_ext(oss.str(), hv_time);
		exts[i]->addExt(new func_const(-0.55));
		exts[i]->addExt(new func_sine(0.5, 2.0, 0));
		exts[i]->addExt(new func_response(-3.0));
	}
	
	for( int i=0; i<nneurons.get(); i++ )
	{
		for( unsigned int j=0; j<nconnections; j++ )
		{
			unsigned int nsrc = lrand48() % nneurons.get();
			real power = drand48() * 7.0 - 4.0;
			real delay = drand48() * 1.25 + 0.25;
			ntime_t psi2 = drand48() * 3.0 + M_LN2 / hv_time;
			ntime_t psi1 = M_LN2 / hv_time;
			tsynapse_messfunc *sy = new tsynapse_messfunc(*exts[i], delay, new func_exp_diff(power, psi1, power, psi2, -1e+30), mess);
//			sy->setDeb(true);
			exts[nsrc]->addSynapse(sy);
		}
		scheduler.scheduleEvent( drand48()*10, makePulse( * exts[ i ], 5.0 ) );
		scheduler.scheduleEvent( drand48()*simulate_until, makePulse( * exts[ i ], 5.0 ) );
//		exts[i]->addSynapse(new tsynapse(obs, 0.1, 1.0));
//		exts[i]->setDeb(true);
	}

	if( nlogging.get() > nneurons.get() )
		nlogging.set(nneurons.get());
	cout << "Logging " << nlogging.get() << " out of " << nneurons.get() << " neurons." << endl;
	for( int i=0; i<nlogging.get(); i++ )
		logger.add(*exts[i], 0.0, false, false, i*4 );
	
	cout << "Produce GNUPLOT plotting file '" << plotfilename << "'." << endl;
	ofstream ofs2(plotfilename);
	logger.gnuplot_def(ofs2, datafilename);

///////////////////////////////////////////////////////////////////////
	cout << endl << "Start Test..." << endl;

	scheduler.run(simulate_until);

	cout << "Test Finished." << endl << endl;
///////////////////////////////////////////////////////////////////////

	cout << "Total fire observed: " << obs.getNPulses() << endl;

	cout << "Total fires: " << totalfire << endl;
	cout << "Total pulses: " << totalpulse << endl;
	cout << "Total partitions: " << totalpartition << endl;
	cout << "Total peak search (peak not found / peak found / non-convex case): " << totalpeaksearch[0] << "/" << totalpeaksearch[1] << "/" << totalpeaksearch[2] << endl;
	cout << "Total peak enclosing: " << totalpeakenclosing << endl;
	cout << "Total partition no-newton (0th/1st/2nd/delta): " << totalpartition_nonewton[0] << "/" << totalpartition_nonewton[1] << "/" << totalpartition_nonewton[2] << "/" << totalpartition_nonewton[3] << endl;
	cout << "Total partition    newton (0th/1st/2nd/delta): " << totalpartition_newton[0] << "/" << totalpartition_newton[1] << "/" << totalpartition_newton[2] << "/" << totalpartition_newton[3] << endl;
	cout << "Total re-scheduled by new pulses: " << totalrescheduled << endl;
	cout << "Total filtered maxgrad:      " << totalfiltered_maxgrad << endl;
	cout << "Total filtered incontinuity: " << totalfiltered_incontinuity << endl;
	cout << "Total filtered nextpulse:    " << totalfiltered_nextpulse << endl;
	return 0;
}