Differences

This shows you the differences between two versions of the page.

--- examples:sim_poisson [2013/12/09 14:50] – created zenke
+++ examples:sim_poisson [2016/02/04 17:39] (current) – typo zenke
@@ Line 1: / Line 1: @@
-====== The important bits ======
+====== sim_poisson example ======
-Since Auryn simulations are executable c programs there is quite some overhead (see below). However since the header logic is usually the same let's focus first on the interesting bits, where the interesting stuff happens.
+Simulates an array of Poisson processes to illustrate the use of [[manual:PoissonGroup]]. The output is written to a [[manual:ras]] file in the same directory.
-<code cpp>
+Full source code: [[https://github.com/fzenke/auryn/blob/master/examples/sim_poisson.cpp]]
-PoissonGroup * poisson = new PoissonGroup(size,kappa);
-poisson->seed(seed);
+===== Running code =====
-sprintf(strbuf, "%s/%s.%d.ras", dir.c_str(), file_prefix.c_str(), world.rank() );
-SpikeMonitor * smon_e = new SpikeMonitor( poisson, strbuf, size);
-sprintf(strbuf, "%s/%s.%d.prate", dir.c_str(), file_prefix.c_str(), world.rank() );
-PopulationRateMonitor * pmon_e = new PopulationRateMonitor( poisson, strbuf, 1.0 );
-RateChecker * chk = new RateChecker( poisson , -1 , 20.*kappa , 10);
-if (!sys->run(simtime,false))
-	errcode = 1;
-</code>
-Running this script should usually in a very brief output like this one
+Running the example program ''sim_poisson'' yields an output like this one
 <code shell>
-bash-4.1$ ./sim_poisson --seed 1231231
+zenke@foobar:~/auryn/build/examples$ ./sim_poisson
-seed set to 1231231.
+[=========================] 100%      t=1.0s  f=5.1 Hz
-[=========================] 100%      t=10.0 s  f=21.3 Hz
 ( 0) Freeing ...
 </code>
-In addition to this the script will create a [[manual:ras]] and a [[manual:prate]] output file (''poisson.0.ras'', ''poisson.0.prate'') as well as a log file ''poisson.0.log''.
+In addition the script will create a [[manual:ras]] and a [[manual:prate]] output file (''poisson.0.ras'', ''poisson.0.prate'') as well as a log file ''poisson.0.log''.
-====== The full script ======
+===== Visualizing the output =====
+Our simulation has written its output (1 second of Poisson spiking from 1000 Poisson neurons firing at 5Hz) to ''poisson.0.ras'', a human-readable [[manual:ras]] file which contains the Poisson spikes. You can peek into the file using a text editor. However, to visualize the  data you need a plotter. I generally like using [[gnuplot]], but any other plotting software which allows reading time series data from columnar ASCII files will do (e.g. [[matplotlib]] in conjunction with [[numpy]]).
+If you have gnuplot installed you can take a quick peek using the command line
+<code>
+echo "plot 'poisson.0.ras' with dots lc rgb 'black'" | gnuplot -p
+</code>
+from within gnuplot
+<code gnuplot>
+plot 'poisson.0.ras' with dots lc rgb 'black'
+</code>
+Either way the output will look similar to this:
+{{ :examples:poisson_output.png |}}
-<code cpp>
-#include <iostream>
-#include <iomanip>
-#include <stdlib.h>
-#include <string>
-#include <boost/program_options.hpp>
+===== Randomness and seeding the random number generator =====
-#include <boost/mpi/environment.hpp>
-#include <boost/mpi/communicator.hpp>
-#include <boost/mpi.hpp>
-#include "auryn_global.h"
+If you run the code again, you will notice that Auryn will generate the same spikes. This is due to the fact that it uses pseudo random numbers which is good for debugging. If you want to generate a different Poisson spike trains try running the code with a random seed. In this simulation this can be done with a command line parameter, but typically this will be done automatically somewhere in your code.
-#include "auryn_definitions.h"
-#include "System.h"
-#include "Logger.h"
-#include "PoissonGroup.h"
-#include "SpikeMonitor.h"
-#include "PopulationRateMonitor.h"
-#include "RateChecker.h"
-using namespace std;
+<code shell>
+zenke@foobar:~/auryn/build/examples$ ./sim_poisson --seed 182319
+[=========================] 100%      t=1.0s  f=5.1 Hz
+( 0) Freeing ...
+</code>
-namespace po = boost::program_options;
-namespace mpi = boost::mpi;
-int main(int ac, char* av[])
+===== Some more options =====
-{
-	string dir = "./";
+This simple simulation has some more parameters. To see which ones they are, you can call the script with
-	string file_prefix = "poisson";
+<code shell>
+zenke@foobar:~/auryn/build/examples$ ./sim_poisson --help
+Allowed options:
+  --help                produce help message
+  --simtime arg         simulation time
+  --kappa arg           poisson group rate
+  --size arg            poisson group size
+  --seed arg            random seed
+</code>
+As you can see you can easily change a few key parameters of the simulation such as runtime, or the Poisson firing rate. It's always a good idea to export some paramters of your code through command line arguments. How this is done will become clearer in the examples. Things that require less flexibility, such as the structure of your simulation itself will be fixed in the simulation file, which is a file with C++ code. Let's have a look at this code in the present example now.
-	char strbuf [255];
+====== The code behind all this ======
-	string msg;
-	NeuronID size = 1000;
+Let's look at the code inside.
-	NeuronID seed = 1;
+You find the source for simulations in the directory ''examples/sim_poisson.cpp''. Since Auryn simulations are executable C++ programs there is quite some overhead (see below). However since the header logic is usually the same let's focus first on the interesting bits, where the interesting stuff happens. The interesting stuff starts here:
-	double kappa = 5.;
-	double simtime = 10.;
-	int errcode = 0;
+<code cpp>
+PoissonGroup * poisson = new PoissonGroup(size,kappa);
+poisson->seed(seed);
+</code>
+The the program first creates a pointer of type ''PoissonGroup'' and assigns an instance of [[manual:PoissonGroup]] to it. The parameters that are passed with the constructor are the size and the firing rate of the group.
+The call of the member ''seed'' seeds the random number generator of all PoissonGroups in a simulation, although in the present case there is only one. With these two calls our 1000 Poisson neurons are ready to be simulated and generate Poisson distributed spikes at the rate ''kappa'' (5Hz in the current simulation). To _see_ these spikes we have to either add other network components or a monitor that writes the spikes to file. We do the latter by issuing the following commands:
-    try {
+<code cpp>
+sprintf(strbuf, "%s/%s.%d.ras", dir.c_str(), file_prefix.c_str(), world.rank() );
+SpikeMonitor * smon_e = new SpikeMonitor( poisson, strbuf, size);
+</code>
+which creates a [[manual:SpikeMonitor]]. The ''sprintf'' command is  used to put together the output filename and store it in ''strbuf''. To construct the actual [[manual:SpikeMonitor]] we simply have to provide the source from where to record (for us that is the poisson group which we find via its pointer ''poisson''), the filename where to write the spikes to (stored in ''strbuf'') and an optional parameter specifying from how many neurons to record. Here we use ''size'' which basically means we are going to record from all neurons. In a next step we will create a [[manual::PopulationRateMonitor]]:
-        po::options_description desc("Allowed options");
+<code cpp>
-        desc.add_options()
+sprintf(strbuf, "%s/%s.%d.prate", dir.c_str(), file_prefix.c_str(), world.rank() );
-            ("help", "produce help message")
+PopulationRateMonitor * pmon_e = new PopulationRateMonitor( poisson, strbuf, 1.0 );
-            ("simtime", po::value<double>(), "simulation time")
+</code>
-            ("kappa", po::value<double>(), "poisson group rate")
+Here the logic is similar as before. The rate monitor needs a source and a filename. The third parameter in this case specifies the binning of the population rate in seconds and it is optional.
-            ("size", po::value<int>(), "poisson group size")
-            ("seed", po::value<int>(), "random seed")
-        ;
-        po::variables_map vm;
-        po::store(po::parse_command_line(ac, av, desc), vm);
-        po::notify(vm);
-        if (vm.count("help")) {
-            cout << desc << "\n";
-            return 1;
-        }
-        if (vm.count("kappa")) {
-            cout << "kappa set to "
-                 << vm["kappa"].as<double>() << ".\n";
-			kappa = vm["kappa"].as<double>();
-        }
-        if (vm.count("simtime")) {
-            cout << "simtime set to "
-                 << vm["simtime"].as<double>() << ".\n";
-			simtime = vm["simtime"].as<double>();
-        }
-        if (vm.count("size")) {
+Before we finally get to run the simulation we have to specify a [[manual:checker]]. Checkers are special lightweight online-monitors. They can communicate with the [[manual:System]] simulation environment and interrupt a run if something goes horribly wrong. This secenario is not unlikely when putting plasticity in recurrent neural networks. Although strictly speaking we do not have to expect this from our simple Poisson group, it is good practice to already get used to it here. We define a [[manual:RateChecker]] which will kill a run if the rate of ''poisson'' rises above 20*kappa or drops below -1 (since this is obviously not possible for a firing rate this a cheap way of sayign "I do not want a lower threshold"). Last but not least the 10 tells the checker to average the rate over a 10s time window. This is also precisely the reason why when running the simulation the value displayed in the end on the command line ( f=21.3 Hz ) is much higher than the actual firing rate of 5Hz. With one second run time the internal moving average did not have enough time to reach the actual mean value.
-            cout << "size set to "
+<code cpp>
-                 << vm["size"].as<int>() << ".\n";
+RateChecker * chk = new RateChecker( poisson , -1 , 20.*kappa , 10);
-			size = vm["size"].as<int>();
-        }
-        if (vm.count("seed")) {
-            cout << "seed set to "
-                 << vm["seed"].as<int>() << ".\n";
-			seed = vm["seed"].as<int>();
-        }
-    }
-    catch(exception& e) {
-        cerr << "error: " << e.what() << "\n";
-        return 1;
-    }
-    catch(...) {
-        cerr << "Exception of unknown type!\n";
-    }
-	// BEGIN Global stuff
-	mpi::environment env(ac, av);
-	mpi::communicator world;
-	communicator = &world;
-	sprintf(strbuf, "%s/%s.%d.log", dir.c_str(), file_prefix.c_str(), world.rank());
-	string logfile = strbuf;
-	logger = new Logger(logfile,world.rank(),PROGRESS,EVERYTHING);
-	sys = new System(&world);
-	// END Global stuff
-	PoissonGroup * poisson = new PoissonGroup(size,kappa);
-	poisson->seed(seed);
-	sprintf(strbuf, "%s/%s.%d.ras", dir.c_str(), file_prefix.c_str(), world.rank() );
-	SpikeMonitor * smon_e = new SpikeMonitor( poisson, strbuf, size);
-	sprintf(strbuf, "%s/%s.%d.prate", dir.c_str(), file_prefix.c_str(), world.rank() );
-	PopulationRateMonitor * pmon_e = new PopulationRateMonitor( poisson, strbuf, 1.0 );
-	RateChecker * chk = new RateChecker( poisson , -1 , 20.*kappa , 10);
-	if (!sys->run(simtime,false))
-			errcode = 1;
-	logger->msg("Freeing ...",PROGRESS,true);
-	delete sys;
-	if (errcode)
-		env.abort(errcode);
-	return errcode;
-}
 </code>
+Finally the simulation is run by issuing ''run'' the following way
+<code cpp>
+if (!sys->run(simtime,false))
+	errcode = 1;
+</code>
+here ''run'' takes at least one parameter which is the simulation time ''simtime''. That is one second in this example. The second parameter is optional. The ''false'' in the example tells run to ignore the RateChecker that we painfully added in the last step. This does not make much sense here, but neither did the Checker in the first place. We will see that the switch is useful in other examples such as [[sim_background]], where one can avoid checking during a transient priming period.
+The full source code can be found here [[https://github.com/fzenke/auryn/blob/master/examples/sim_poisson.cpp]]